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We present the theoretical and analytical bases of optimal techniques to measure 
weak gravitational shear from images of galaxies. We first characterize the geometric 
space of shears and ellipticity, then use this geometric interpretation to analyse images. 
The steps of this analysis include: measurement of object shapes on images, combining 
t~-L ■ measurements of a given galaxy on different images, estimating the underlying shear 

from an ensemble of galaxy shapes, and compensating for the systematic effects of im- 
age distortion, bias from PSF asymmetries, and "dilution" of the signal by the seeing. 
These methods minimize the ellipticity measurement noise, provide calculable shear 
uncertainty estimates, and allow removal of systematic contamination by PSF effects 
to arbitrary precision. Galaxy images and PSFs are expressed as "Laguerre expan- 
sions," a 2d generalization of the Edgeworth expansion, making the PSF correction and 
shape measurement relatively straightforward and computationally efficient. We also 
discuss sources of noise-induced bias in weak lensing measurements — selection biases, 
and "centroid" biases arising from noise rectification — and provide a solution for these 
H . and previously identified biases. 

Subject headings: gravitational lensing — methods: data analysis — techniques: image 
processing 



1. Introduction 

Gravitational lensing is a powerful tool for studying the distribution of matter in the Universe, 
because photons are deflected by all forms of matter, regardless of luminosity, composition or dy- 
namical state. Dramatic manifestations of lensing — multiple images, Einstein rings, and giant arcs, 
so-called strong lensing — provide much information on the highest overdensities in the Universe, 
namely rich galaxy clusters, cores of individual galaxies, or collapsed objects. To characterize the 
more typical mass structures, or those without a fortuitously aligned bright background source, we 



may use weak gravitational lensing, in which we analyze the low-order distortions of the ubiquitous 
background galaxies in order to infer the mass distribution. Weak gravitational lensing signals 
are extraordinarily subtle, even by astronomical standards: one seeks a shear (or magnification) 
of the galaxy images amounting to a few percent at most, more typically 0.2-1% in current stud- 
ies. Because the undistorted image is not observable, the lensing distortions must be detected as 
a perturbation to the intrinsic distribution in galaxy shapes (or sizes), which have variation of 
30% or more, giving a signal-to-noise ratio (S/N) of ~ 1/30 from observation of a single galaxy. 
Hence a very large number of galaxies must be observed before the weak lensing becomes detectable 
over this intrinsic shape noise. Weak lensing analyses could not even be attempted until automated 
means of measuring very large numbers of galaxy shapes became available (Valdes, Tyson, & Jarvis 
1983; Tyson et al. 1984). Furthermore, optical and atmospheric distortions in a typical sky image 
cause coherent shape (and size) distortions that can masquerade as a lensing signal. Such system- 
atic errors are 1-10% in a typical image, up to 50 times larger than the weak lensing signals. A 
means to remove this contamination is crucial; the necessary analyses can only be conducted with 
well-calibrated, linear detectors. 

Successful detection of a weak lensing signal did not occur until CCD images of sufficient depth 
and field were available (Tyson, Valdes, & Wenk 1990), and early detections were of the rs 10% 
shears that are found in the inner regions of rich clusters of galaxies (Fahlman et al. 1994; Bonnet, 
Mellier, & Fort 1994; Smail et al. 1995). In regions of strong shear, the S/N is sufficiently high 
that a map of the lensing mass can be created (Kaiser, Squires, & Broadhurst 1995). Mellier 
(1999) includes a review of results from cluster-lensing studies. 

With the increase in collecting area of CCD imagers, sufficient background galaxies can be 
measured to allow convincing detection of smaller shear signals around weaker overdensities: around 
individual weak clusters (Fischer et al. 1997) or collections of galaxy groups (Hoekstra et al. 2001); 
around individual galaxies (Fischer et al. 2000; Smith et al. 2000; Wilson et al. 2001b). Most 
dramatically, lensing signals on random lines of sight, caused by the background matter fluctuation 
spectrum, have now been detected and are one tool for "precision cosmology" (Wittman et al. 2000; 
Van Waerbeke et al. 2000; Bacon et al. 2000; Wilson et al. 2001a). As technology has advanced, 
weaker and weaker shears have become detectable under the shape noise, sometimes as small as 
a few tenths of a percent (e.g. Fischer et al. 2000, Jarvis et al., in prep.). As a consequence, 
the demands for rejection of systematic errors have become more stringent. In many current weak 
lensing publications, it is clear that the uncorrected systematic effects are only slightly smaller than 
the signals under study. It is therefore fair to say that, at present, it is the analysis techniques, 
rather than the ability to collect galaxy images, that bar the way to higher precision in many 
weak-lensing studies. 

This paper describes the techniques for extraction of weak-lensing signals from imaging data, 
which we have developed over the past few years to meet these increasing demands. As described 
below, our efforts focus on the shear rather than the magnification of the galaxy images by the lens, 
and hence we are measuring galaxy ellipticities. The desiderata for a weak-lensing methodology 



include: 

1. Shapes of individual galaxy images are determined with the highest possible accuracy in the 
presence of measurement noise on the image. 

2. Each shape measurement should have a known error distribution. 

3. Individual galaxy shapes should be combined to yield an estimate of the underlying lens shear 
with maximal S/N. 

4. The shear estimator should have an error level and a calibration that can be derived directly 
from the data, without recourse to Monte-Carlo simulations. 

5. The galaxy shapes should be corrected for the systematic biases due to the point-spread 
function (PSF), to arbitrary precision. 

6. The scheme must allow for a PSF that varies continuously across the image and is different 
in each exposure. 

Given the intrinsic floor on weak lensing accuracy because of shape noise, one might ask why 
we should expend much effort on goal (1), which is to minimize the effects of measurement noise — 
normally, we consider that once the ellipticity measurement noise a e is <C <jsn = (e 2 /2) 1 ' 2 ~ 0.3, 
further gains do not increase the shear estimation accuracy — the error erg on the lensing distortion 
8 will just become (tsn/vN, with N the number of measured galaxies. We note first that the 
sky density of galaxies scales with apparent magnitude m as 10 am with 0.3 < a < 0.4. If we 
can cut the shape measurement error for a given image noise level, then we can either use fainter 
galaxies in our lensing measurement (increasing N), or cut the required exposure time. Second, 
note that convolution with a PSF suppresses the measured lensing signal and the intrinsic shape 
noise. Hence the level to which we aim to reduce a e must, for poorly resolved galaxies, be well 
below the canonical 0.3. Thirdly, we will see in §4.2 that it may be possible to measure the shear to 
an accuracy much better than asN /\N , in cases where the distribution of intrinsic galaxy shapes 
in the ellipticity plane has a cusp or pole at e = 0. In simple cases such as a population of circular 
disk galaxies, the accuracy to which we can measure the applied shear can increase without limit 
as the measurement noise is decreased. 

The need for traceable uncertainties is also critical as weak lensing is used to measure the power 
spectrum of mass fluctuations in the Universe. In this application, the measurement uncertainties 
(including shape noise) contribute to the power spectrum and must be accurately estimated and 
subtracted to reveal the true cosmic power spectrum. Of course an accurate calibration is also 
necessary for most applications to precision cosmology; if one must rely on simulated data for the 
calibration, there is always the danger that the simulations do not properly incorporate some aspect 
of the real world. 



Finally, the need for removal of the systematic PSF ellipticities to arbitrary precision is ex- 
tremely strong. In the course of this paper we will try to describe other approaches to the problem 
and compare to our own. The methods in most common current use (e.g. Kaiser, Squires, & Broad- 
hurst (1995), KSB) are formally valid only certain special cases of PSF. While heuristic adjustment 
and testing has demonstrated that the method works to nominal accuracy in more general cases, 
the absence of a generally valid method is troubling. A formally exact PSF correction scheme has 
been put forward by Kaiser (2000) [K00], which is based upon a Fourier-domain calculation of the 
effects of shear and of PSF convolution. Our approach will be to decompose the image and the 
PSF into a vector over orthogonal polynomials, and treat the deconvolution as a matrix operation 
carried to desired order. A very similar approach has been independently put forth by Refregier 
(2001). 

This is a longwinded paper, likely to be read in detail only by practitioners of weak lensing. A 
more casual reading will be beneficial to those who wish to understand the methods and limitations 
of past and future weak lensing analyses. Some of the techniques we develop may be useful beyond 
the weak lensing analysis, for example our deconvolution method (§6.3.5) and the methods for rapid 
convolution with spatially varying kernels (§7). As discussed by Refregier (2001), our orthogonal- 
function decomposition can be a useful means for compression of galaxy images. 

The paper outline is as follows: the following section describes the mathematical space oc- 
cupied by ellipticities and shears. Understanding the geometry of this space makes it easier to 
see how our (and other) measurement schemes work. In §3, we describe a scheme which uses our 
geometrical conceptualization of ellipticity to produce measurements with maximal S/N in images 
with infinitesimal PSF; a formula for the resultant uncertainty in each ellipticity is also derived. 
Next, §4 discusses several schemes for combining shape measurements of a given galaxy from differ- 
ent exposures and/or filter bands, to obtain the shape estimate that again offers the best possible 
S/N and a closed-form error estimate. §5 describes the means to combine shape estimates from 
different galaxies to form an optimal estimate of the underlying lensing shear. In the absence of 
measurement noise, this takes a simple closed form; in the presence of measurement noise, some 
approximations must be made to obtain a closed form for the calibration and error of the shear, 
and hence we do not fully satisfy goal (4) above. §6 is a very extensive discussion of the effects of 
the PSF on the image, other approaches to the problem, and our method for optimal extraction of 
the intrinsic shape in this case. In this section we introduce the Laguerre decomposition technique. 
§7 uses the Laguerre formalism to construct convolution kernels that can add symmetries to the 
PSF of an image; this is one means of removing the ellipticity biases due to the PSF. §8 discusses 
two very important effects that can give rise to biased lensing measurements even when a perfect 
deconvolution for PSF effects is possible. It is likely that these biases are present in all previously 
published data. 

Finally, §9 puts together all of the methods developed in the paper in a flowchart form describ- 
ing how raw image data are converted into optimized, calibrated lensing shear data. We reserve for 
a succeeding paper (Jarvis et al. 2002) the detailed discussion of the code that implements these 



methods, and a verification of its performance on real and simulated data. In Appendices to this 
paper, we present the formulae for invoking various transformations on the Laguerre-decomposition 
representation of an image, and derive some approximate PSF-correction formulae that were used 
for the analyses of Smith et al. (2000) but which are superseded by the full Laguerre methodology. 



2. Geometry of Shape and Shear 

2.1. Linear Approximation to Lensing 

The goal of weak gravitational lensing studies is to infer a distant gravitational potential via 
the distortions that the potential's deflection of light imparts upon the population of galaxies in the 
background. The lensing is fully characterized by the map u(x) from the observed angular position 
x to the source angular position u. The surface brightness observed at x is equal to that which 
would have been observed at u in the absence of the lens. For an individual background galaxy 
that is not near a lensing caustic, the map can be accurately approximated by a Taylor expansion 

u = — x + u . (2-1) 

ax 

The displacement uo carries no information (unless the source is multiply imaged) because the 

source plane is unobservable. The 2x2 amplification matrix has a unique decomposition of the 

form 

rfx „„ 

— = ^SR, 2-2 

du 

where R is an orthogonal matrix (rotation) ; S is a symmetric matrix with unit determinant (shear) ; 
and \x is a scalar magnification. 

The rotation R is not useful for lensing studies because the unlensed orientation is not known, 
and the ensemble of background galaxies should be isotropic and hence any collective statistic 
should be unchanged by rotation. Furthermore the rotation is absent in the limits of single-screen 
or weak lensing. 

The magnification \x increases the angular size by /j, and the galaxy flux by fi 2 . While the 
unlensed quantities are not observable, the magnification is still detectable because the mean flux 
and size of the population will shift. The magnification also reduces the sky-plane density of sources 
by fi 2 . The magnification thus modulates the number vs magnitude relations for a given class of 
background galaxies, in a manner which depends upon the size/magnitude/redshift distribution of 
the original population. 

The shear S has two degrees of freedom, corresponding to the ellipticity and position angle 
imparted on a circular source galaxy. For weak lensing this shear is undetectable on a single 
galaxy because the unlensed shape is not necessarily circular and is not observed. The collective 
distribution of galaxy shapes is assumed to be intrinsically isotropic, and the applied shear breaks 
this symmetry, rendering it detectable and measurable. 



Both the shear and the magnification thus produce measurable effects on the ensemble of 
galaxies and can in theory be used to quantify the potential. Shear measurements have been 
used for numerous quantitative studies, but magnification methods still yield at best marginal 
detections (Dye et al. 2001). There are several factors that favor the shear method: first, the two 
effects of magnification (increased flux and reduced areal density) push the counts of background 
sources in opposite directions, weakening the signal. More importantly, the shear is manifested as a 
variation in the mean orientation of galaxy shapes, and this mean is zero in the absence of lensing; 
the magnification signal is a modulation of N{m) or some other non-zero quantity. It is always 
far easier to measure a small change from zero than a small change in a non-zero quantity. For 
example, exploitation of the magnification effect in the weak-lensing regime would require absolute 
photometry to much better than 0.01 mag accuracy. Magnification measurements, on the other 
hand, give a direct measure of the projected mass, whereas mass reconstructions from shear data 
are degenerate under the addition of a constant-density mass sheet. Hence magnification data are 
very useful when there is no a priori means of determining the mean mass overdensity in the image. 

Henceforth we will ignore the magnification effect and describe how to optimally measure the 
shear S. 



2.2. Parameterizations of Shear 

2.2.1. Diagonal Shears 
The simplest shear matrix is a small perturbation aligned with the coordinate axes: 

The effect of this transformation upon a circular source-plane object is to induce an elongation 
along the x axis, creating an elliptical image with axis ratio q = b/a = 1 — r\. We can use this 
matrix as a generator for the full family of diagonal shear matrices with arbitrary rj to obtain 




S v =[ /2 , -oo < n < oo. (2-4) 



The set of diagonal shear matrices forms a group under simple matrix multiplication. The operation 
is commutative, and clearly corresponds to simple addition of the r) parameters: 

s m = s ri2 x s m <=^ V3 = m + m- (2-5) 

For this reason we will call r) the conformal shear and will find it a useful parameterization of shear. 
Other common parameterizations of shear include the axis ratio q, the distortion 5 (Miralda-Escude 
1991), and the reduced shear g = 7/(1 — k) (Schneider & Seitz 1995), which are related to r/ via 

q = b/a = e~\ (2-6) 



a 2 



a 2 + b 2 



t&nhr], (2-7) 



9 = ^=tanhr//2. (2-8) 

Bonnet & Mellier (1995) define a further set of shear parameterizations, also easily expressed in 
terms of r/: 

e = \-q = 2e~' ?/2 sinh7//2 
I -a 2 



l + q 2 
l + q 2 



tanh?7 



cosh?? 



2q 
t = e5b = sinh 77. 

Note that for rj <C 1, the parameters rj, e,£,r, and the distortion <5 are all equal. Note also that 
most other author's formulations of shear do not define the matrix S to have unit determinant, so 
do not form a group. 



2.2.2. General Shear 
A general (non-diagonal) shear matrix can be decomposed into a diagonal shear and rotations 



as 



« _n CD __ | COSn 2 + Sin ^ 2 COS ^ Sin ^ 2 Sm ^ 1 fO Q\ 

S>„ )/9 - K^K^ - ^ g . nh 2 s . n fl cogh 3 _ ginh 2 cog J • i^-9) 

S^ transforms a circular source to an ellipse with axis ratio q = e~ v at positional angle (3 = 9/2. 
The shear can be represented as a 2-dimensional vector 

r] = (rj + ,r) x ) = (rj cos 9,rj sin 6). (2-10) 

Likewise a shear may be represented as a two-dimensional distortion (<5 + ,<5 x ), etc. A shear (r/ + ,0) 
creates ellipses oriented to the x or y axes, while (0, rj x ) aligns circular sources to axes at 45° to 
the coordinate axes. The shear rj is not a vector in the image space, but rather is a vector in a 
non-Euclidean shear space that we describe below. 

The full set of shear matrices do not form a group under matrix multiplication as S-n S-n may 
be asymmetric (two-screen lenses can effect a rotation for this reason). But we can form a group 
with an addition operation for 2-dimensional shears defined as 

*73 = *72©*7i «=► S *?3 R = S *7 2 S *7i' (2" 11 ) 

where R is the unique rotation matrix that allows S/i to be symmetric. The geometric meaning 
of the shears is preserved since R will leave a circular source unchanged. The simplest expression 



of the composition operation in terms of components is 

cosh 7/3 = cosh r/2 cosh r/i + sinhr/2 sinhr/i cos (^2 — ^l) (2-12) 

sinhr/3 sin(03 — $2) = sinh?7i sin(0i — 82). 

Note that the second equation is not symmetric in the two operands and hence the shear matrix 
group is non-Abelian. The identity element is r] = and the inverse of t] = (r/+,rj x ) is —t] = 
(— r/ + , — f] x ). The addition formula in terms of distortion components is derivable from (2-12), and 
is given by Miralda-Escude (1991): 



^3+ = ; , x x > ( 2 " 13 ) 



<Ji+ + <W + (W^)[l - VI - S'i](S lx d2+ - 


- s 1+ s 2x ) 


1 + <$i • s 2 




Six + S 2x + (S2+/5'i)[l - v 7 ! - Wi+^x - 


- <5ix<52+) 



03X : 1 + <W 2 

We omit the derivations of these equations, which are straightforwardly but tediously executed by 
composing the transformation matrices. A more elegant derivation follows from noting that the 
transformation Equation (2-9) transforms the complex plane as 1 

z ^ cosh-z + e ie sinh-z. (2-14) 

It will be useful to consider the limit where 62 <C 1: 

(dd©o) x « dd x j 

If we instead make 5\ <C 1, the asymmetry of shear addition is manifested as a change to the 
azimuthal component formula: 

(*ed«) x « VT^d6 x tW<1 ' 5x=a (2 " 16) 



2.3. The Shear Manifold 

We define a metric distance between two points 773 and rj 1 in shear space as 

s=\V3~Vi\ = \m\, V2 = V3®{~Vi)- (2-17) 

The differential form of the metric can be derived by specializing Equations (2-12) to the case 
ds = r]2 < 1, 0i = 0, 2 = 0, yielding 

m = 771 + ds cos (2-18) 

3 tanh??! = (is sin (2-19) 



We thank the anonymous referee for this derivation. 



which means that the metric is 

ds 2 = ( m -r jl ) 2 + tmh 2 m {e 3 -e 1 ) 2 (2-20) 

= drj 2 + tanh 2 77 d9 2 (2-21) 

= {l-S 2 ) 2 dS 2 + S 2 d6 2 . (2-22) 

Note that the r\ version of the metric has the normal Euclidean form for the radial component, 
and the 5 parameterization has the Euclidean form for the tangential component of the metric, but 
neither representation gives a fully Euclidean metric — the shear space is curved. The 2-dimensional 
shear manifold defined by this metric can be embedded in Euclidean 3-space as illustrated in 
Figure 1. This geometric depiction of shear is helpful in understanding the transformations of 
shears. Near r/ = the surface is tangent to the Euclidean plane, so small shears add with Euclidean 
component- wise addition. The shear-space surface then curves upwards and as the conformal radius 
rj grows large, the surface approaches a cylinder of unit radius about the z axis. If we project the 
shear surface onto the z = plane, the radius vector in this plane is equal to the distortion 5. The 
S vector is confined within the unit circle. 



2.4. Definition of Shape 

A shear is a transformation of the image plane; we next need a quantity to describe the shape 
of an arbitrary galaxy image. Let / represent some object whose isophotes are a family of similar 
ellipses. We can simply parameterize the shape of / by the shear rj which produces this object from 
some object 1$ having circular isophotes, i.e. 

I-q = Sfjlo =>• S r)Ir) = Iq- (2-23) 

We could thus call rj the conformal shape of the object, and can think of a given ellipse as a location 
on the shear manifold. More commonly the distortion is used to define the shape; an object is said 
to have ellipticity e if a shear with distortion 8 = — e makes it circular. We will use the symbol 
e since this quantity agrees with the traditional second-moment definition of ellipticities for truly 
elliptical objects. Equation (2-23) makes it obvious how an ellipse with shape r] 1 will be transformed 
under the action of a shear rj 2 : we simply add the shear to the shape using the addition rules of 
shear space [Equation (2-12)]: 

Sj ?2 Ir 7i = /? ?3> with ^3 = V2 © Vi- ( 2 " 24 ) 

Likewise we can also say that a distortion S maps the ellipticity e — ► S © e. In general, an applied 
shear may be viewed as a shift of all shapes along the shear manifold. We will use e to represent 
the shape of an image, whereas S represents a shear, which is a transformation of the image plane. 
The shape and shear spaces, however, transform identically under an applied shear. 

A real galaxy has some image intensity distribution /(x) which may not have elliptical isophotes; 
we would like to define a shape for an arbitrary image. By analogy to Equation (2-23), all that is 
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Fig. 1. — The shaded surface is an embedding of the shear manifold into Euclidean space. The radius vector 
along this surface is the conformal shear rj; the radius upon projection onto the xy plane is the distortion 
5 (or ellipticity e). At small r\ the manifold is tangent to the S plane, and at large 77 approaches the unit 
cylinder. Two shear vectors si and S2 of length r\ = 1 are plotted from the origin, both on the true shear 
manifold and in the S plane. The result of adding the two vectors is also plotted; displacements do not 
commute in this non-Euclidean space. 
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needed is some definition of a "round" image. Let M(J) be any measurement applied to the image 
which has the simple property that if M(7) = 0, then M(SrtI) 7^ for any non-zero shear 77. Then 
the condition M(J) = is our definition of roundness, and we can assign a unique shape rj to an 
arbitrary object I(x) by the condition 

M(S_ T? I) = 0. (2-25) 

Any shape defined by such a rule clearly transforms under an applied shear just as an ellipse does, 
namely via Equation (2-24). With any roundness measure M we therefore have a definition of 
shapes and their mapping under shears that follow the rules of addition in shear space. We do not 
attempt to prove that the solution to Equation (2-25) exists or is unique. 



2.5. Shear in Fourier Space 

For consideration of the effects of convolution upon sheared images it will be useful to ponder 
the action of shears in Fourier space. We first note that shearing an image I(x) by rj is equivalent 
to shearing its Fourier transform i"(k) by —77. For a diagonal shear: 

S^I(k x ,k y ) = (27T)- 1 fd 2 xl(e- r >/ 2 x,e r >/ 2 y)exp[2m(k x x + k y y)} (2-26) 

= (27T)- 1 f d 2 xI(x,y)exp[2TTi(e^ 2 k x x + e^ 2 k y y)] (2-27) 

= I(e^ 2 k x ,e-^ 2 k y ). (2-28) 

A nondiagonal shear must also satisfy this relation since rotation of the real-space function corre- 
sponds to rotation in /c-space. We can therefore just as easily define a galaxy's shape by a roundness 
criterion in fc-space as in real space. This is useful when considering finite resolution (§6) or when 
analyzing interferometric images with limited Fourier coverage. 



3. Optimal Measurements (without Seeing) 

In this section we derive an optimum method of measuring galaxy shapes in the case where 
the angular resolution and sampling of the instrument are assumed to be perfect. In §6 we will 
treat the more complex case of seeing-convolved images. 



3.1. The Ideal Test for Roundness 

We have defined the shape of an image in Equation (2-25) by asking what coordinate shear 
is needed to make the object appear round. We need to choose a measurement M which detects 
with the highest possible signal-to-noise ratio (S/N) any small departure of the image from its 
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round state M = 0. It can be shown that, under some sensible simplifying assumptions about 
M, the solution of Equation (2-25) becomes equivalent to finding the best least-squares fit of an 
clliptical-isophote model to the galaxy image. 

We assume first that the measurement will be a linear function of /(x). Any non-linear method 
will prove extremely difficult to apply to the case where the image has been convolved with a PSF— 
it will be hard to use the measurements of bright stars to correct the shapes of faint galaxies for 
convolution. The most general form of M is then 



M(/) = /d 2 xW(x)/(x), (3-1) 



where W is some weight function. The weight will be two dimensional, as the measurement must 
test for departure from roundness in both the rj + and r/ x directions in shear space. 

We consider first the weight component to detect a small change in r] + . We can decompose 
the image I(r,9) into multipole elements I m (r) via 

oo 

I(x)=I(r,e) = Yl ^(rY™ 9 (3-2) 

m=— oo 

I m (r) = ±- [ n d9I(r,e)e- me (3-3) 

^ Jo 



We are interested in the change in our measurement upon mapping of the image /(x) — > /(x) = 
/(Sx) where S is a shear of amplitude r] <C 1 oriented on the x-axis, as in Equation (2-3). The 
quantity for which we wish to optimize the S/N can thus be written as 

00 />oo 

SM = M{I) - M(I) = Y, / rdrw m {r)[I m {r)-I m {r)], (3-4) 



m=— oo 



with w m an arbitrary radial function for each multipole. A little bit of algebra yields the transfor- 
mation of multipoles 



L(r) - I m (r) = - A [(m- 2)/ m _ 2 - rl' m _ 2 - (m + 2)/ m+2 - rl' m+2 ] + 0( V 2 ), (3-5) 



V 

4 

where the primes denote derivatives with respect to r. For an object with truly circular isophotes, 
we have I m {r) = for m ^ 0, and the only effect of the shear is to induce a quadrupole term 
I 2 = —r]rlQ/4. For objects without perfect circular symmetry, there are terms beyond the monopole. 
But for an object to be thought of as "round," the monopole term Iq should dominate the higher 
multipoles. The monopole is also the only term guaranteed to be positive for all galaxies. Hence the 
largest effect of the shear will be to alter the m = ±2 quadrupole intensities (which are conjugates 
of each other as / is real) . The optimal sensitivity to small shear should therefore weight only the 
quadrupole term: 

M(I) = J r dr I 2 (r)r 2 w(r) = — f / r dr d8I(r,9)w{r)r 2 e~ 2ie '. (3-6) 
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It is clear that this quadrupole test is the optimal linear measurement for objects with circular 
symmetry; for more general shapes, the shear has effects on other multipoles that can be measured 
and used to enhance signal-to-noise (related to the suggestions of Refregier (2001)). But this would 
require a knowledge of Ia{v) and the other multipoles to construct the ideal formula; we settle on 
the simple quadrupole as the best general solution, as we are always guaranteed that I is present 
and positive for any real source. The measurement of some weighted quadrupole is also the normal 
definition of ellipticity for weak lensing measurements (Miralda-Escude 1991). 



Combining Equations (3-4), (3-5), and (3-6) we obtain 

-n 



5M 



__]_ J r drr w(r) [rl' + 4/4 + rl'4] . 



(3-7) 



The noise in the measure of M can be derived in two limits: the most common case will be sky- 
dominated observations, for which the variance of the flux in area A is nA, where n is the number 
of sky photons per unit area (it is assumed that / is in units of photons). In this case we have, 
from Equation (3-6), the variance of each component of SM 



Var(M) 



n 
4^2 



r dr dOr w (r) cos 29. 



(3-8) 



If we ignore the I4 terms in Equation (3-7) as being dominated by Iq terms, then the choice of 
weight function which optimizes the detectability of the shear n is 

ldlo 

r dr 



ir, 



opt I 



< 



(3-9) 



With this optimal weight, the variance of the measurement M would lead to an error in each 
component of r/ equal to 
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(3-10) 



If the object is much brighter than the night sky, then the noise is no longer uniform and the 
optimization becomes 

4^y j 

Id In Jo 



Var(M) 



rdrdOr w (r)I(r) cos 26 



(3-11) 
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3.1.1. Gaussian Objects 



An elliptical Gaussian object, when sheared to be circular, will obey 
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(3-12) 
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In the sky-limited case the optimal weight is the same Gaussian: 

w opt (r) = ^o,e- r2 / 2 °\ (3-13) 

r 

Note that the optimal weight for shape measurement is in this case equal to the optimal filter for 

detection, i.e. a matched filter. If we define the detection significance v as the signal-to-noise ratio 

for detection of the object with the matched filter, we find 

^ = [JdAw(r)I(r)] 2 
JdAnw 2 (r) 

f 2 

1 (3-15) 



Amna 2 



2 lGrrna 2 

°r l = — 7T~ ( 3 ~ 16 ) 



2 x2 



(3-17) 

We therefore end up with the simple result that the error in each component of the shear is 2 over 
the detection significance. 

The above derivations assumed that the center and the size a of the Gaussian were known in 
advance. If there were a sky filled with Gaussian galaxies, we likely would not know in advance 
the size and location of each. We can determine the centroid in the usual manner by forcing the 
weighted first moments to vanish: 

dAwIre w = 0. (3-18) 

The weight for centroiding does not necessarily have to match that used for the shape measurement, 
but it is convenient to do so. The proper size a w for the weight can be forced to match the size a 
of the object by requiring the significance to be maximized: 



dv 



[dAwI(l-r 2 /a%). (3-19) 



In the limit of a Gaussian with low background noise, the Equations (3-11) apply and the 
optimal weight is uniform. The detection significance in this case is just v = \fj (with the flux 
/ in photons), and we find again that the standard error in rj is equal to 2/ v. In practice this 
situation can never be realized because the weight extends to infinity, and at large radii the sky or 
read noise will dominate the shot noise from the galaxy, and neighboring objects will impinge upon 
the integrations. We will henceforth confine our discussion to background-limited observations. 



3.1.2. Exponential or Other Profiles 

To obtain an optimized measurement of a real galaxy we would have to measure its radial 
profile and construct a custom optimized weight using Equation (3-9). The majority of galaxies 
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are spirals or dwarfs which are typically described by exponential profiles: 

I(r)o,e~ ar =* Wopt ( r ) = ^-. (3-20) 

This weight diverges at the origin, though all the necessary integrals of the weight are convergent. 
If the galaxy is truly cusped in the center, then the intensity near the center is very sensitive to 
small shears and is weighted heavily. 

In practice it is simpler to adopt a weight function that is universal (up to a scale factor), 
especially at low S/N where attempts to measure each individual profile would be pointless. There 
are a number of reasons to prefer a Gaussian weight: 

• The Gaussian drops very quickly at large radii, minimizing interference from neighboring 
objects. Integrals of all moments are convergent. 

• Weights with central divergences or cusps are difficult to use in data with finite sampling, 
and also amplify the effects of seeing on the galaxy shapes. The Gaussian is flat at r = 0. 

• Gaussian weights are analytically convenient, allowing many useful formulae to be rendered 
in closed form. 

• Gaussian weights allow construction of the family of orthonormal basis functions that we will 
use in later chapters to compensate measured shapes for finite resolution. 

• The Gaussian is not far from optimal for most galaxy shapes. For a well-resolved galaxy 
with an exponential profile, the Gaussian weight measures r] with only 7% higher noise than 
the optimal weight in Equation (3-20). In the presence of seeing, the difference between the 
Gaussian and optimal weight is even smaller. Recall that any weight we choose yields a valid 
definition of roundness and hence of shape; the Gaussian just incurs a small penalty in noise 
level. 

The procedure for measuring galaxy shapes is therefore as follows: 

1. Estimate a shape r\ for the image I and apply the shear S_« to obtain /. 

2. Iterate the center and size of the Gaussian weight function until the centroid condition (3-18) 
and the maximum-significance condition (3-19) are satisfied. 

3. Compute the second moments with the Gaussian weight function 

M(I)= [dAw(r)I(r,0)r 2 e- 2ie , w(r) = e'^^. (3-21) 



4. If the real and imaginary parts of M are zero, then ij is the shape of the object. If not, then 
we use the measured M to generate another guess for rj and return to step 1. 
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The process is mathematically equivalent to measuring the second moments of /(x) with an el- 
liptical Gaussian weight, and iterating the weight ellipticity, center, and size until they match the 
measured object shape. It is therefore an adaptive second-moment measurement. The method is 
also mathematically equivalent to finding the elliptical Gaussian that provides the best least-squares 
fit to the image. 



3.2. Uncertainties in Shape Estimates 

Once we have settled upon a weight of the form w = e~ r ' 2u , we can integrate SM from 
Equation (3-7) by parts and use (3-8) to calculate the variance in rj. We will again ignore the I& 
terms; this means we may have a small tendency to over- or under-estimate our shape errors if 
galaxies tend to be boxy or disky. We first define the weighted flux f w , significance v, and weighted 
radial moments {Ir m ) w as 



f w = dAwI (3-22) 

v 2 = /2/Var(/ w ) = fl/nna 2 (3-23) 

(Ir m ) w = -!- [dAwIr m . (3-24) 

Jw J 

The condition (3-19) for optimal significance is 

(Ir 2 ) w = a 2 (3-25) 

and under this condition the variance in each component of the shape is 

* = S& +0( ° (3 - 26) 

A 



i/ 2 (l -a 4 ) 2 



+ 0(0 (3-27) 



a 4 ee ^-1. (3-28) 

The quantity 04 is a form of kurtosis which is zero for a Gaussian image. The terms of order i/~ 4 
arise from errors in the centroid determination, and are discussed further in §8. 

The procedure for measuring an object of shape rj requires applying a shear —rj to the image 
coordinates x to produce a coordinate system x' in which the object appears round. The uncer- 
tainties in Equation (3-26) apply in this sheared coordinate system. Because the object is round in 
this frame, there is no preferred direction in the shear space and the uncertainty region is circular, 
with an uncertainty of ±<r^ on each component (r]' + ,r]' x ) measured in the sheared frame. We must 
reapply a shear +77 to restore our measurement to the original coordinate system. This process 
is illustrated in Figure 2. The 77 coordinate transformations defined in Equations (2-12), in the 
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limit of r/i = ffq < 1, indicate that the uncertainty region on r\ will be elliptical, with a shrunken 
principal axis in the circumferential direction of the shear manifold: 

oq = — ^ — ^> ds = t&nhr] (Tq = <r„sechry. (3-29) 

smhr/ 

If we instead use Equations (2-16), we can find the uncertainty ellipse in the ellipticity plane to be 

ae = ( l^> (3-30) 

eae = v 1 — e <t»j. 

So on the e unit circle, the uncertainty ellipse shrinks radially by 1 — e 2 and tangentially by yl — e 2 
as we transport the error region from the origin back to the original ellipticity e. 

Note that our derivation assumes that the noise characteristics of the image are unchanged 
when we apply a shear. This is true in the background-limited case because the noise spectrum 
is flat, and our shear matrices S have unit determinant. For an image that has been smoothed or 
deconvolved, the power spectrum will have structure and the shear will alter the noise statistics. 
We will discuss this in §6 in the context of finite image resolution. 



3.3. Comparison with Other Methods 

Galaxy ellipticities for weak lensing were first determined by computing unweighted second mo- 
ments of the intensity (Tyson, Valdes, & Wenk 1990). If the moment integrals are taken to infinity, 
then the measured ellipticities transform under shear using the addition rules in Equation (2-13), 
and furthermore the correction for PSF effects is extremely simple. It is clearly impractical, how- 
ever, to carry the integrals to infinity, since neighboring objects will interfere, the noise is divergent, 
and, as noted by K00, many common PSFs have divergent second moments. So the initial methods 
generally used some sort of isophotal cutoff to the moments. This has the disadvantage of creating 
moments which are non-linear in the object flux. An alternative would be to use unweighted mo- 
ments within a fixed circular aperture, but as noted by KSB, the noise properties of unweighted 
moments are far from optimal. 

In the KSB method, the measured ellipticity exsB is computed from the second moments 
measured with a circular Gaussian weight with size selected to maximize the detection significance 
v. [A different weight function was suggested by Bonnet & Mellier (1995).] The distinction between 
the KSB method and ours is that KSB always apply a weight which is circular in the original image 
plane; in our adaptive method, the weight is circular in the sheared image plane which makes the 
object round. Or, as viewed in the original image coordinates, the weight is an ellipse with shape 
iterated to match that of the object. This distinction has two consequences: first, our adaptive 
method yields lower uncertainties for non-circular objects because the weight is a better match to 
the image. This effect is minor, though, for objects with e < 0.5, and a minority of images arc 
more elliptical then this. 
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Fig. 2. — The scheme for ellipticity measurement and its errors are illustrated schematically: in (a), the 
triangle marks the true shape of a target galaxy, located in the e plane. The shape is determined by shearing 
the image until the galaxy appears round. In the ellipticity plane, we are moving the object along the vector 
to the origin. Panel (b) shows the location of the target (and the original coordinate grid) in the e plane 
after application of the shear that makes it round. The shaded region represents the uncertainty region for 
the shape in the sheared view — the error region must be circular because the image is round. Finally in (c) 
we undo the applied shear, shifting the target and the error region back to the original shape. This mapping, 

3 2 ) in the radial (azimuthal) axes. 



however, shrinks the error region by a factor 1 



(vT 
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The second, more important advantage over the KSB method is that our definition of shape via 
Equation (2-25) guarantees that our measured ellipticities e are transformed by an applied shear 
via Equations (2-13). The circular- weight exsB does not have this property — indeed for an object 
with elliptical isophotes, exsB does not equal the true ellipticity. In the KSB method, it is therefore 
necessary to calculate a "shear polarizability" for each object, describing the response of exsB to a 
small shear 5; this polarizability depends upon the radial profile and m = 4 moments of the object. 
The "polarizability" of our measured ellipticities is just the i5<1 limit of Equations (2-13): 



e + d ^e + S-e(S-e). (3-31) 



1 + 6-t 

This transformation rule, and the shapes of our uncertainty regions in ellipticity space, arise solely 
from the geometry of the shear manifold and are independent of the details of the galaxy images. 
This will simplify the following discussions of methods to derive a shear from an ensemble of 
measured galaxy shapes. 

We implement the adaptive-weighted-moments scheme in the program ellipto, described 
further in Jarvis et al. (2002). 



4. Combining Exposures 

In a typical observing program, a given background galaxy is imaged in a number of different 
exposures, in one or more bandpasses. This is done to increase exposure time, permit rejection 
of cosmic rays, and/or gather color information. Multiple exposures can also reduce systematic 
effects by placing data for a given galaxy on different parts of the detector and in different seeing 
conditions. 

We hence encounter the question of how to combine data on a given galaxy in different images 
to an optimal single measure of the shape. There are two possible approaches: 

1. Measure the shape on each exposure, then create a weighted average of the measurements as 
the final shape. 

2. Register and average the images, then measure the object on the combined image. 

We first consider which offers the lowest noise on the final shape. Consider the task of combining 
iV exposures, with the object having significance uq on each exposure. Following Equation (3-26), 
the uncertainty in the shape of a nearly-round object measured from a single image will be 

<o = 4 + C" [Var(x )] 2 = 4u 2 + CV ~ 4 . (4-1) 

The second term is the uncertainty due to centroiding error, and C and C' are constants of order 
unity. If we average measurements (Method 1), then we decrease a v by y/N. If we average images 
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(Method 2), then we increase v by \[N . The net error on the shape in the two cases is then 



_ y /1 + C/2ul Method 1 (42) 

Nu I y/1 + C/2Nv$ Method 2 

The two methods are equivalent, except for the centroiding noise. If uq is not > 5 than averaging 
images will produce better accuracy on ij. Keep in mind that (a) the galaxy will not even be 
detected on the individual exposures unless vq > 3, and (b) the galaxy is useless for weak lensing 
unless a v < 0.5, which requires Nvq > 16. When N < 3, the centroiding penalty is small for 
any object that will be useful, so a combined image is extraneous. When N > 5, there are many 
galaxies detectable on the summed image that are not detectable on the individual images, and a 
summed image has detectability and centroid-noise advantages. 

There is a compromise, "Method 1.5," which has the practical advantages (delineated below) 
of Method 1, while retaining the small S/N edge of Method 2: that is to create a summed image 
and use it for object detection and centroid determination, so that Var(xo) ~ 1/Ni/q. Then this 
centroid is used to measure shapes on individual exposures, and the shape errors are equivalent 
to Method 2. In practice we will combine deconvolved Laguerre coefficients (§6.3) rather than 
measured shapes. 

There are several reasons why it may be preferable to average catalogs instead of images: 

• Correction of shapes for PSF effects is paramount, and only possible if the PSF is constant 
or slowly varying across the image. If the different exposures in a summed image overlap 
only partially, then the PSF (and noise level) will jump discontinuously as one crosses the 
boundaries of component exposures. It is therefore preferable to correct for PSF effects on 
an exposure- by-exposure basis. If the PSF is very stable (e.g. a space telescope) or if the 
exposures all have nearly the same pointing (a single deep field), then a summed image will 
have well-behaved PSF variations. 

• For the smallest objects, the exposures with the best seeing will contain nearly all the useful 
shape information and should be weighted heavily (Kochanski & Fisher 1994). Large objects 
are, on the other hand, measured equally well in every exposure. Averaging catalogs allows 
one to adjust the weights of different exposures on an object-by-object basis, whereas this is 
not possible when combining images. 

• If there are exposures in different bands, then the optimal weighting of the exposures is 
dependent upon the color of the object. This is easily done when averaging cataloged shapes 
but not easily done by summing images. 

• Creation of a summed image requires registration and interpolation of pixels. The latter 
process smoothes the noise field and causes subtle variations in the PSF, both of which 
complicate later analyses. 
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• An especially pernicious hazard to creating a summed image is that slight mis-registration of 
the component images will cause coherent elongations of the images, which if not corrected 
will mimic a lensing signal. This is discussed by KSB; in theory such effects are handled by 
a proper PSF correction scheme. This is a danger for Method 1.5 as well. 

Some practical advantages to Method 2, combining images, are that: 

• The data storage and processing requirements can be lower for a single combined image if iV 
is large. 

• In Method 1, outliers (from cosmic rays) are rejected on an object-by-object basis, whereas in 
Method 2 rejection is pixel- by-pixel. If the galaxies are very oversampled and the cosmic-ray 
rate is high, Method 2 could salvage the un-contaminated parts of galaxy images that Method 
1 discards. 

For the simplest circumstances (a single-filter stack of images with common pointing), image av- 
eraging is easier and has few drawbacks. For multi-filter or mosaicked data, catalog averaging is 
needed. The hybrid Method 1.5 is best for such cases, though more work. In the rest of this section 
we detail procedures for each Method. 



4.1. Combining Images 

There are standard tools for combining exposures into a single image. We remark here upon 
a few special considerations when doing this for weak lensing observations. 

First, accurate registration is paramount. Our scheme for image registration is described in 
Jarvis et al. (2002). 

Second, the use of median algorithms is commonplace but dangerous. Proper correction for 
PSF effects will require that the images of bright stars have precisely the same PSF as do the faint 
galaxies. But with a median algorithm, the bright, high-S/N stars will be constructed with a PSF 
which is a median of all the exposures. The images of faint objects, however, will tend toward 
a PSF that is the mean of all the exposures, because the noise fluctuations will dominate PSF 
variations. The final PSF will therefore vary with magnitude. A sigma-clipping average is much 
preferred over the median for the necessary task of cosmic-ray rejection when combining. 

Similarly one must be careful about rejecting saturated pixels. There will be many stars which 
saturate only on the best-seeing exposures; if the saturated pixels are rejected, these stars will have 
final PSFs which are broader than the PSF for faint objects. One must take care to ignore stars 
which are saturated in any one of the exposures. 
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4.2. Combining Shape Measurements 

Suppose that a given galaxy has been measured to have ellipticity e^ in images i € {1,2,..., N}. 
We desire the e which best estimates the true ellipticity of the object. Using the results of §3.2, we 
see that in the absence of PSF distortions, the minimum-variance estimate of e will be that which 
minimizes the \ 2 given by 

x 2 = ^ (eee . )S -l (eee .) (4 . 3) 

i 

Here, e e^ is equivalent to e © (— e^), where © corresponds to the addition operator introduced 
in Equation (2-11), and Sj is a covariance matrix which is cr 2 I in simple cases. 

Note that if the e^ are measured in different filters, than the galaxy may have no single well- 
defined ellipticity. By "best estimate," then, we must mean that which offers the best sensitivity 
to a weak lensing distortion, and the minimum-variance combination of the ej is still the desired 
quantity. 

The is a non-linear operator, so we could use a non-linear minimization algorithm to find 
the value of e at which \ 2 ls minimized. However, this is both impractical for time considerations 
and unnecessary since the values of a v are usually small. Thus, we can linearize the subtraction 
operator 

e©e i = T i (e-e i ) + 0((e-e i ) 2 ) (4-4) 



T can be derived from Equations (2-13) 
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(4-5) 



The linearized \ 2 becomes 



which has a minimum at 
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(4-6) 
(4-7) 

(4-8) 



This is a standard least-squares solution for the mean of the ej given covariances Xlj in a 
Euclidean e space. In the simple case of S = a^I, the expression for Sj simplifies considerably to 

1 



^(1-e 2 ) 2 



1-e 2 e x e + 
e x e+ 1-e 2 , 



(4-9) 
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which is equivalent to Equation (3-30). However, when we apply corrections for PSF dilution, we 
will find that the covariance matrix is more generally an ellipse with axes aligned in the radial and 
tangential directions. That is, 

£ = R, ( °l 2 ° 2 ) R e (4-10) 

where R# is a rotation matrix with 9 = arctan f — J . In either case the minimization (4-8) is 
numerically straightforward, and we are left with an uncertainty ellipse for the mean e. It is wise 
to implement some outlier-rejection algorithm in this process as well. 



5. Estimating Shear from a Population of Shapes 

Now we presume to have measured ellipticities e^ for a set of TV distinct galaxies, with known 
measurement uncertainties for each. Our final task is to create a statistic 8 from the ej which best 
estimates the lensing distortion 8 that has been applied to this ensemble. There are three main 
effects which must be considered in constructing the estimator: first, the e^ respond differently to an 
applied distortion 8, as embodied by Equation (3-31) for true ellipticities, or by shear polarizabilities 
for the KSB estimators, so we need to know the responsivity 1Z = 88/ 88 of our statistic. Second, 
the variety of ellipticities in the parent (unlensed) galaxy population causes shape noise in the 
shear estimate. In most weak lensing projects this is the dominant random error, and we wish to 
minimize its effects. Third, there is measurement error in each ellipticity, which we also wish to 
minimize in our shear estimator. 

Most practitioners have adopted a simple arithmetic mean of e + and e x as estimators for 
the applied distortion {e.g. Fischer & Tyson 1997). Using the weak-distortion Equations (3-31), 
it is easy to see that, in the absence of measurement error, this estimator has a responsivity 
1Z = 1 — <7g N and a variance Var(<5 + ) = a^ N /N, where we have defined the shape noise <t| n = (efj.) 
(the x component has the same properties and shape noise). 

Others have realized, however, that rare, highly elliptical galaxies have too much influence on 
the arithmetic mean and should be deweighted. Cutoffs on |e| (Bonnet & Mellier 1995) or other 
weighting functions w(e) (Van Waerbeke et al. 2000) have been applied to the ellipticities and 
tested with simulations, but without any sort of analytic optimization or justification. Lombardi &: 
Bertin (1998) consider the optimization of a general weighted sum of second moments (rather than 
ellipticities); this unfortunately couples the ellipticity measurement to the distribution of sizes of 
the galaxies and leads them to consider only weights which are power-law functions of the moments. 

Hoekstra et al. (2000) present a weighting scheme which incorporates both measurement error 
and the shape noise, and K00 gives a detailed discussion of optimal weighting for distortion mea- 
surements. Both are similar to our method in many respects, which we comment upon at the end 
of this section. 
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Fig. 3. — The left panel shows a model for intrinsic distribution P(e) of galaxy shapes over the e plane — it 
must have circular symmetry. When each galaxy is sheared by 5+ — 0.05, the galaxy distribution shifts to 
the right as in the middle panel. The right-hand panel shows the change in population under the applied 
distortion; this is the signal which we wish to detect. Shape noise arises from the Poisson fluctuations in the 
population, which is proportional to the left-hand panel's -P(e). The optimal weight for <5 + determination is 
the ratio of the right-hand to the left-hand panel. 

5.1. Without Measurement Error 



We start with an unlensed background galaxy population with ellipticities distributed within 
the unit circle according to 

.7 V 

P(e')d 2 e' = P(e')e'de'de'. (5-1) 



dN 



A fundamental assumption of weak lensing is that the background is isotropic so that the unlensed 
population can have P depend only upon the amplitude of e', not its orientation. The effect of 
a distortion 6 is to map the background population to a new, anisotropic distribution P$(e), as 
illustrated in Figure 3. We are given a sample of iV galaxies from the new distribution, and our 
task is to estimate the 6 which gave rise to the distribution from the original P. 

One approach is to find the value of S which maximizes the likelihood of the observed ej. This 
is true when 

° = ^E lo g^^) = M;E lo g p (i-^© e ^- ( 5 - 2 ) 

A similar condition holds for <5 X • These equations define the maximum-likelihood S even for strong 
distortions — though there is not in general a closed- form solution for strong S. 

For weak lensing (<5 <C 1), Equation (3-31) and the conservation of number can be used to 
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derive Ps(e) to first order in 5. 
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It can further be shown that the maximum-likelihood estimator for 8 takes the form 

l-e 2 dlogP' 



-y< 



dc 



(5-3) 
(5-4) 

(5-5) 

(5-6) 

(5-7) 

(5-8) 



The parenthesized expression is thus a weight function for combining ellipticities into a distortion. 
We can also show that this weight function is optimal in terms of S/N for weak distortions, as 
follows. Let us create an estimator S which is a general weighted sum of the ellipticities, 



s= E w(ej)ej = Jd 2 ew(e)eP 5 (e) 
J2 w ( e i) fd 2 ew(e)P(e) 

The response of this statistic to a small applied shear is 



(5-9) 
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(5-10) 



(5-11) 



where in the last line, we have dropped terms linear in e+ or e x which average to zero over an 
isotropic population. With an isotropic population, the derivative <9<5 X /dS x = 1Z as well, and the 
off-diagonal elements of dS/dd are zero. 

We may use Equation (5-11) to calculate the response of any weighted estimator by summing 
over the observed e,, because the small difference between observed and intrinsic distributions does 
not alter 1Z to first order. In the case where we have some analytic form for P(e), we may replace 
the sums with integrals over the distribution to obtain 
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(5-12) 
(5-13) 
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In the absence of measurement noise, the variance in S is due to shot noise. Assuming that 
the background galaxies obey Poisson statistics and their shapes are randomly assigned, we can 
propagate the Poisson errors through Equation (5-9) to get the expected error 
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fd 2 ew 2 (e)e 2 + P(e) 
N[Jd 2 ew(e)P(e)}' 
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(5-14) 
(5-15) 



In (5-15) it is assumed that the sum is over a sufficiently large ensemble of background galaxies 
to sample the distribution -P(e). Any weak lensing measurement has thousands of background 
galaxies, so this gives a direct estimate of the error in the shear. 



The optimal weight is that which minimizes Vai(5+)/lZ 2 , which is 
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(5-16) 

(5-17) 

(5-18) 



where the last line gives the optimized error in 5+/1Z, which is our calibrated estimate of the 
distortion. Equation (5-16) reproduces the maximum- likelihood solution in Equation (5-8). This 
may be compared to the distortion uncertainty for equal weighting w = 1, 
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We first see that a simple arithmetic mean of the ellipticities is the optimum estimator only if 
P oc (1 — e 2 ) a for some exponent a. For the real galaxy population, there can be a significant gain 
in accuracy through use of w op t over equal weighting. An extreme case is a population of randomly 
oriented circular disks, for which 
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(5-20) 
(5-21) 



With w = 1, we would have as = 0.590/-^^. The optimal weight diverges at e — ► to take 
advantage of the extreme sensitivity of -P(e) to distortion near e = 0. The integral in Equation (5- 
18) in fact diverges at e = 0, driving as to zero — which would be a significant improvement over the 
equal-weighting case! Unfortunately any small measurement error or departures from circularity 
for the disks will smooth out the central spike in -P(e), creating a finite value for a$- 
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Figure 4 shows the P (e) measured for well-measured galaxies in the CTIO lensing survey 
(Jarvis et al, in preparation). These shape histograms are derived from 230,000 galaxies which 
are well resolved (R > 0.4, under the definitions in Appendix C) and have errors on the intrinsic 
ellipticity of o r] < 0.03 — primarily galaxies of magnitude 17 < jur < 22. The shape of P{e) is 
observed to be highly dependent upon the surface brightness (SB) of the galaxies. The low-SB 
galaxies show the rise at e > 0.8 expected of a disk population, but there the distribution drops at 
e > 0.95 instead of diverging — this reflects the finite thickness of the disks. There is also no pole 
at e = for the low- SB galaxies, showing that the disks are not perfectly circular. The high- SB 
galaxies are presumably early types since there are very few with e > 0.5. While the value of 
P(e = 0) increases with surface brightness, it always remains finite, but with dP/de < 0. The ideal 
weight Equation (5-16) therefore grows as 1/e as e — ► 0, but the contribution to the S/N does not 
diverge at zero as for perfect disks. None of the Pie) curves is well fit by a single Gaussian or 
power law. 

The intrinsic ellipticity variance <tsn varies from 0.30 to 0.49 between the highest and lowest 
SB bins. The optimally weighted distortion S/N per galaxy for high-SB, early types is 2-3 times 
higher than for the lowest-SB galaxies, indicating the desirability of incorporating some galaxy- 
type discriminant — surface brightness, color, or concentration — into the weighting scheme. The 
only requirement is of course that the discriminant be independent of ellipticity. It seems likely 
that -P(e) will vary substantially with magnitude. 

The heavy histogram in Figure 4 combines all well- measured galaxies with 20 < m^ < 21, 
which we henceforth use as a representative measure of the real galaxy population. The distribution 
has ctsn = 0.38, which would lead to as = 0.44/yiV for an unweighted average. The optimal 
weighting gives ag = 0.33/viV; the weighting therefore gives a gain equivalent to a 1.8-fold increase 
in N. The gain in telescope time is at least as large. This gain is reduced, however, in the presence 
of measurement noise, which will tend to wash out the sharp feature in P(e), as discussed next. 

We reiterate two favorable results of this section: first, the responsivity 1Z and the variance of 
8 can be expressed exactly as direct sums over the observed population, for arbitrary choice of w; 
there is no need for calculation of polarizabilities nor recourse to simulated images. Second, we note 
that the variance in 5 can be significantly below the canonical crg N /N if the ellipticity distribution 
P(e) has structure that is not washed out by measurement noise. 



5.2. With Measurement Error 

A galaxy image with true ellipticity e will be measured at some ellipticity e, with a probability 
distribution of p(e|e). We consider a population of galaxies all having the same significance v and 
resolution parameter R (see §6 and Appendix C), so that they all have a common p(e|e). The 
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Fig. 4. — The distribution of intrinsic ellipticities for modestly bright galaxies (rriR rs 20) is plotted; we 
plot 2ireP(e) rather than P(e) as the nature of the population is more apparent. The distribution is highly 
dependent upon surface brightness fj,R, presumably reflecting the difference between spheroid- and disk- 
dominated galaxies. The dashed line is the distribution for an isotropic population of circular disks. The 
high-/i galaxies are more useful for distortion measurements. The heavy histogram combines all surface 
brightnesses in the magnitude range 20 < wir < 21. Though it is difficult to tell from this plot, P(e) is 
finite and increasing as e — > 0. Optimal weighting takes advantage of this structure to reduce the noise in 
the distortion measurement. 
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measured distribution of ellipticities observed under distortion S will then be 

P S (e) = Jd 2 eP s (e)p(e\e), (5-22) 

where Px is the distribution of true ellipticities as in the previous section. The symmetry of P and 
p in the ellipticity plane guarantees that in the absence of distortion, the measured distribution 
P{e) must again depend only upon the magnitude, not the direction, of the measured ellipticity. 

Equation (5-22) is not strictly a convolution because the measurement function p(e\e) may 
depend upon e and not simply upon (e — e) — for example Equation (3-30) describes how the error 
ellipses contract as e departs from zero, even if the significance of the detection is held fixed. In §6 we 
will show that the behavior of the measurement error is different when the effects of PSF smearing 
upon the image are important, so at this point we will consider p(e|e) to be, most generally, some 
kind of Gaussian whose 1<7 ellipse depends only upon the magnitude e. An important point is that 
the functional form of p(e|e) is unchanged by an applied distortion, since p is determined by v and 
R, which are unchanged by a pure shear. 

Another fact to keep in mind is that, with finite resolution and noise, it is possible to measure 
e > 1, if the image noise makes the object appear smaller than the PSF in some dimension. 
Our formulae should therefore be tractable even for e > 1, and we cannot simply discard such 
measurements without contemplating the consequences. 

We proceed as in the previous section, by assuming a distortion estimator of the form 

~ _ IXeQej = fd 2 ew(e)P s (e)e 

£«>(ei) J<Pew(e)P(e) ' 

(n) = dS+ fd 2 e [w(e)e+ fd 2 ep(e\e)(dP s (e)/d6 + )] 

[ ' d5+ jd?€w{e)P{e) ( ' ' 



fd 2 e w{e)e + fd 2 eP{e)p(e\e)e+ (3 



l-e 2 dlogP 
e de 



JcPew(e)P(e) 



(5-25) 



Var( . +) _ i^mHs, (5 _ 26) 



N[f<Pew(e)P(e)] 2 
J2w 2 {e)e 2 + 
E>(e)] 2 



(5-27) 



™°pt( e ~) = ^~ fd 2 eP(e)p(e\e)^ U- 1 ^ dl J P ) (5-28) 



P(e) J e + \ e de 

Given functional forms for the intrinsic distribution P(e) and the uncertainty function p, we could 
use (5-28) to derive an optimal weight, and use (5-25) and (5-26) to get the responsivity and noise 
for the estimator using this or any weight function. In most cases these integrals will not have 
analytic solutions. 
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The parenthesized quantity under the integral in Equation (5-13) is a galaxy's responsivity to 
shear, which depends upon the intrinsic shape. Equation (5-25) is the average of this responsivity 
for the galaxies with some measured shape. The measurement noise can cause these two quantities to 
differ; in other words, a naive determination of the responsivity is biased by the measurement noise. 
KSB-based methods will also suffer a calibration error due to this effect; binning the polarizabilities 
in parameter space can reduce the noise in the polarizability, but will not remove biases. Precision 
cosmology will require that such calibration issues be addressed — there are unfortunately no cosmic 
calibration standards for shear. 

The need for P(e) in the above formulae is an unfortunate complication since P is the directly 
observed quantity. Note that the variance of the estimator can be expressed as a closed sum over 
the observed shapes (5-27), but the responsivity cannot. A precise calibration of the resultant 
shear/mass maps requires, therefore, that P be estimated either by deconvolving the observed P 
with the error distribution p, or by recourse to higher-quality images that give P directly. 

Derivation of the the optimal weight also requires knowledge of the intrinsic P, but we can 
explore some generic cases, and make some approximations that give workable methods. 



5.2.1. Approximate Form For Responsivity with Errors 

We wish to have a form for 1Z as a sum over the observed objects and applied weights, as 
in Equation (5-11), for the case of finite measurement errors. Toward that end we can take the 
derivative of Equation (5-23), which is greatly simplified if we assume that the measurement error, 
i.e. §j — ej, does not have any first-order dependence on 5. While not strictly valid it is a good 
approximation. In this case 



£ w(l - <4> g ) + \% (1 - (e 2 + ) e - - (e + e x ) e -e x /e H 

n = — 



(5-29) 



where the brackets indicate an average of the true quantity at a given measured value, e.g. 

(e+)e = f d 2 ep(e\e)el (5-30) 

J d 2 ep(e\e)P(e)e 2 + 
/d%p(e|e)P(e) (5 " 31) 

Note that the weight function w may depend upon e directly, but also indirectly through some 
dependence in its covariance matrix XI. If the measurement error function p(e|e) and the intrinsic 
distribution P(e) have circular symmetry, then we must be able to write 

<4)g = fco(e) + fci(e)4 (5-32) 

where k$ and k\ are functions only of the magnitude, not direction, of e. We must also have 
(e+e x )e = fei(e)e+e x . Further manipulation, taking advantage of the isotropy of the parent popu- 



lation, yields 

E 

n = — 
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(5-33) 

z> 

This form for TZ depends only upon the observed quantities and the chosen weight scheme, except 
through the two functions ko and hi, which we will approximate below. The resemblance to 
Equation (5-11) is clear. With this equation and some integration by parts, we may also derive a 
form for the optimal weighting function: 

r~ w~M u , ldk ° . ~ dkl l-k -he 2 dlogP 
w opt [e, E(e)] = Sh + -=— + e— - g ^. (5-34) 



5.2.2. Special Case: Gaussians 

In general the functions fco(e) and k\(e) must be calculated numerically using a presumed 
underlying P(e) for the background population, but analytic solutions are possible in the case of a 
Gaussian P(e) with variance o"g N in each component (the shape noise) and a constant measurement 
error a 2 on each component. We find that both ko and k\ are independent of e: 
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(5-35) 



The quantity / is the fraction of the total ellipticity variance that is attributable to shape noise. 
When the measurement noise is small, / ~ 1, the ideal weight is close to (1 — e 2 )/(o"| N + a 2 ). This 
is quite similar to the weight adopted by Hoekstra et al. (2000). 

For the Gaussian case, dlogP/de is also quite simple, so Equation (5-33) can be used. In our 
surveys to date (Smith et al. 2001; Jarvis et al. in prep.) we have adopted a weight that results 
from optimizing the Gaussian case. 



5.2.3. Practical, Nearly Ideal Approximation 

We obtain an approximation to the correct responsivity 1Z and resultant ideal weight if we adopt 
the constant ko and k\ functions in Equations (5-35) even for non-Gaussian P(e) distributions. The 
shape noise Cg N may be defined as the assumed variance of the underlying e + and e x , and may 
be found by subtracting the measurement noise from the observed (e 2 ). The measurement noise 
a 2 is known for each galaxy using the methods of this paper; since the covariance matrix for e is 
generally anisotropic, some representative scalar must be selected. 

With these guesses for ko and k\ in hand, 1Z may be estimated with a sum over the observed 
galaxies using Equation (5-33) for any chosen weight function. 
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For the real-Universe shape distributions measured in the CTIO survey, we find that the 
following "easy" weight function offers very close to optimal distortion measurements: 

w= [e 2 + (1.5^) 2 ]~ 1/2 , (5-36) 

where a^ is the shape uncertainty that the object would have were it circular [cf. Equation (3-30)]. 

We can check the accuracy of our approximations numerically for chosen P(e) and p(e\e) 
functions. We examine the case when P(e) is that shown in Figure 4 for galaxies with 20 < 
run < 21, and the measurement error follows Equation (3-30). We find that the weight function 
given by Equation (5-34) is in fact very close to optimal for all noise levels, even when the simple 
approximations (5-35) are used for the k functions. The "easy" weight Equation (5-36) also performs 
nearly optimally, so most applications could use this weight and need not attempt to determine 
P(e). 

A more critical question is whether the approximations (5-35) yield a proper estimate of the 
calibration factor 1Z when used with Equation (5-33). Figure 5 shows how the simple 1Z estimator 
compares to the correct value in Equation (5-25) for our choice of underlying distribution and 
the "easy" weights (5-36). The approximate form yields a responsivity correct to better than 5% 
f° r (J v ^ 0.4. It is clear from the Figure that some detailed knowledge of the underlying P(e) 
distribution will be needed in order to calibrate lensing measurements to the one percent level. 

The lower panel of Figure 5 shows the potential advantage of optimal weighting. When the 
measurement error is > 0.2, there is little difference between various weighting schemes. For an 
unweighted distortion estimator, the accuracy levels out as the measurement noise drops below 
a e ~ 0.2. When optimal weights are used, however, the distortion errors continue to drop as the 
measurement error is pushed toward zero — the optimal weights take advantage of the e = cusp 
in the shape distribution. Our "easy" weight scheme recovers nearly all of this potential gain. 

To summarize, a practical method of weighting and calibrating the response in the presence 
of measurement noise is: 

1. Determine the underlying <jg N = (e+) of the the intrinsic distribution. The measurement 
error a„ (or the full covariance matrix) of each galaxy is known from the formulae in previous 
sections. 

2. Approximate the quantities ko and k\ with Equation (5-35). 

3. Choose a weight function, for example the "easy" form Equation (5-36), or preferably the 
optimal form Equation (5-34) which can be derived from the observed distribution -P(e). 

4. The distortion estimator is the sum (5-23). The variance in the estimator is the sum (5-27). 

5. The estimator (and variance) must be scaled by the responsivity which is calculated with 
Equation (5-33). 
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Fig. 5. — The upper panel plots the approximate, simplified closed-sum estimate of the calibration factor 
1Z (Equation [5-33]) relative to the exact form (Equation [5-25]), as a function of the ellipticity measurement 
noise at e = 0. The weight function is the "easy" Equation (5-36) and the factors ko and k\ adopt the simple 
heuristic approximation (5-35). The heavy (black) curve is for the population of 20 < rriR < 21 galaxies, the 
upper (green) curve is for a low-SB sample (20 < (j,r < 21), and the lower (red) curve for a high-SB sample 
(17 < (Mr < 19), for which the intrinsic P(e) distributions are plotted in Figure 4. The simple formulation 
yields a calibration accurate to 5% or better in all cases, but 1% accuracy is difficult to achieve. The lower 
panel shows the uncertainty in the distortion determination when the galaxy shapes are combined with 
optimum weights (solid lines), the "easy" weights (5-36, dotted), and equal weighting (dashed). The line 
weights/colors code the galaxy sample, as above. Note that with optimal or easy weighting, the distortion 
errors continue to shrink even when measurement error is well below the canonical shape noise of 0.3. 
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If more is known about the intrinsic shape distribution, then more accurate functions for ko and 
k\ may be derived and used in the sums. 



5.3. Additional Weighting Considerations 

The optimal weight may depend upon parameters other than the observed ellipticity e. It 
must, for example, depend upon the measurement errors as described above. If intrinsic shape 
distribution P(e) depends upon galaxy type, for example, then it may be advantageous to have 
different weight functions for each type — as long as galaxy type can be determined independent of 
e. Ellipticals may be distinguished from spirals, for example, by a concentration ratio such as the 
622 coefficient described below. 

The expected shear will depend upon the redshift z of the background galaxy, and hence there 
may be a z-dependent weight determined from photometric redshift estimates of the background 
population. More crudely the apparent magnitude may be taken as an indicator of z and used 
in the weight formulation. The use of such additional weights depends upon the problem being 
addressed: see Smith et al. (2000) as an example. 



5.4. Relation to Previous Methods 

The optimal weighting scheme of K00 differs from ours in several respects. It assumes a 
different measure of e which does not follow our geometric transformation relations, so the mean 
polarizability of galaxies must be calculated in some parameter space, e.g. of flux, size, and e. 
Binning or smoothing in some such space is common to most of the KSB-based weighting schemes 
as well (Erben et al. 2001; Hoekstra et al. 2000). The variance of the estimators Si are calculated 
within each bin, then weights are chosen inversely to the variance of each bin to create a minimum- 
variance weighted mean estimator. There are three contributors to the variance from each bin: 

1. The intrinsic ellipticities of galaxies within the bin are drawn at random from the parent 
distribution. If |e| is one parameter of the space, then only the direction 6 is allowed to vary 
within the bin. This is of course the shape noise. 

2. The measured shape of a given galaxy is drawn from the measurement-error distribution. 

3. The number of galaxies within the bin is drawn from a Poisson distribution. If |e| is a param- 
eter, then this Poisson noise includes some elements of the shape noise (1) and measurement 
noise (2). 

All three of these effects are important to optimization; all are included in our formulation and 
(implicitly) in that of K00, so we expect them to be essentially equivalent in the long run — this 



35 



is not the case, though, for some of the heuristic or parameter-space weight formulations in the 
literature. The virtue of our scheme is that the nature of the weight function is apparent given 
the intrinsic shape distribution P(e) and the measurement errors, and there is no need to choose a 
parameter space for weight selection and polar izability smoothing. Our formulation tells us when 
further parameters might be desirable, namely when P(e) changes significantly. 



6. Measurements with Finite Resolution 

The preceding sections outline a method for optimal recovery of weak distortion from galaxy 
images, and rigorous estimation of the uncertainties on these shears, for the case when the detector 
views the galaxies with perfect resolution. Unfortunately, the finite resolution of real observations 
has a strong effect upon shape measurements in every weak lensing observation to date, even those 
using the Hubble Space Telescope. Finite resolution produces two deleterious effects: 

1. A PSF which is not circularly symmetric can induce ellipticities on the images, thus breaking 
the intrinsic isotropy of the background galaxy population and mimicking a lensing distortion. 
This is a bias induced by asymmetric PSFs. Since present-day weak lensing surveys are seeking 
distortion signals well below 1%, measured shapes must be compensated for even the slightest 
asymmetry in the PSF with some smear correction. 

2. Convolution by a circularly symmetric PSF will make most galaxies appear rounder, driving 
e — > 0. This is therefore a dilution of the true lensing signal. While this mechanism cannot 
create a lensing signal where there is none, it can misleadingly modulate the lensing signal or 
cause a calibration error in the inferred mass distribution. 

Most (but not all) approaches to PSF corrections treat the bias and dilution effects in separate 
steps. To our knowledge, all published weak-lensing observations have incomplete PSF correction, 
leaving systematic distortion errors of a fraction of a percent or higher. While in most cases these 
residual systematic effects have not altered the validity of the authors' conclusions, the proper 
correction of PSF effects is presently the largest barrier to the use of weak lensing for precision 
cosmology. 

In §6.1 we review some existing approaches to these problems, in §6.2 we contemplate how 
one would ideally wish to approach the problem, and in further sections we develop two means of 
implementing a nearly-ideal approach: one which treats bias and dilution separately, and another 
which corrects both problems simultaneously with a limited form of deconvolution. 
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6.1. Existing Approaches to PSF Corrections 

6.1.1. The Unweighted Ideal 

Unweighted second moments of galaxies are ideal measures of ellipticity — not only do the 
derived ellipticities transform according to the rules of shear space, but correction for PSF effects 
is in principle quite simple because the unweighted second moment of the image is just the sum 
of the original moment and the PSF moment. Thus by subtracting the PSF moments from the 
measured moments we simultaneously correct for bias and dilution, and obtain the image shape. 
In the case where the PSF is round, the dilution of a true (pre-seeing) image-plane ellipticity e* to 
an observed (post-seeing) ellipticity e° is described by the exact equation 

e* = e°/R (6-1) 

R ^ ^ 

(r 2 )i + (r 2 )* 

( r 2\ 

(r 2 )o 

The resolution parameter R is determined by the unweighted second radial moment of the measured 
image (r 2 ) relative to that of the PSF (r 2 )*. Two things to note: first, the error ellipse on the 
dilution-corrected, pre-seeing ellipticity e* is magnified by 1/R from the original measurement 
error Equation (3-30), and is further stretched in the radial direction by the uncertainty in R itself. 
Thus the error ellipse is no longer simply described by a single a v . Second, Equation (6-1) can 
give rise to |e*| > 1, if the noise makes the galaxy look smaller than the PSF about some axis. We 
cannot arbitrarily discard such measurements without creating a bias in our mean shear. These 
two phenomena are common to all modes of PSF-dilution correction. 

This blissfully simple dilution correction is spoiled by two major problems: first, as discussed 
above, unweighted second moments have divergent noise properties and for this and other rea- 
sons are not practical shape estimators. An equally serious problem noted by K00 is that the 
second moments themselves are divergent for many realistic PSFs. Further, many galaxies follow 
deVaucouleurs profiles, for which the second moment converges very slowly. 

The simple formulae (6-1) are still valid under the special circumstances that the object and 
PSF are both Gaussians. The post-PSF object is again a Gaussian, and deconvolution of Gaussians 
is a simple subtraction of second moments. Hence any shape-measuring algorithm which extracts 
the proper ellipticity for a Gaussian ellipsoid would allow PSF dilution correction via Equations (6- 
1) in this limited (and unrealistic) case. 

Some early weak lensing measurements (Valdes, Tyson, &: Jarvis 1983) adopt second-moment 
subtraction as a means of PSF correction, despite the fact that this method is not exact when 
isophotally-bounded or weighted moments are used, and the images are not Gaussian. This would 
not suffice, however, for the more sensitive measurements being done today. 
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6.1.2. Heuristic Methods 

In the case of unweighted moments or Gaussians, Equation (6-1) would indicate that a regres- 
sion of the lensing signal against (r 2 )" 1 would yield a distortion free of PSF effects as (r 2 )^ 1 — ► 0. 
Mould et al. (1994) have attempted to measure very weak shears in the presence of PSF effects 
using such a regression (though against (r)^ 1 , in which case a linear relation is not expected). Even 
with weighted moments, we expect the PSF dilution and bias to decrease as the object becomes 
well-resolved, so there is some basis to this method, even if it is not exact. Other problems, how- 
ever, are that the distortion is not likely to be the same for all sizes of galaxy as they likely lie at 
different distance. Also the regression will lead to substantially higher noise than a more direct 
dilution correction. 

Another approach to the dilution correction is exemplified by FT97, who attempt no analytic 
correction, instead calibrating the dilution effect by measuring simulated background galaxies which 
have been subjected to the same distortion, seeing, sampling and noise as the real images. Such a 
simulation is an essential test of any weak-lensing methodology. The difficulty with sole reliance 
upon simulated data is that the result is extremely sensitive to one's ability to exactly match 
the size-magnitude distribution of the true galaxy population, because the dilution correction is 
a strong function of size (as in Equation [6-1]) in the typical regime of slightly resolved galaxies. 
Further, as we show below, the dilution correction depends upon detailed higher-order moments 
of the galaxy images, which would be very difficult to simulate faithfully. One alternative is to 
use high-resolution, high-S/N images from HST instead of simulated galaxies — but the total sky 
area imaged to sufficient S/N by HST is a tiny fraction of a square degree, too small for rigorous 
calibration tests. It would be preferable to have an analytic correction for dilution, and use the 
simulated data to spot-check the accuracy of the analytic method. 



6.1.3. Perturbative Methods 

A step beyond the unweighted-mean approximation to the bias correction is taken by KSB 
and by FT97. Both make the assumption that the anisotropy of the PSF can be described as 
a small anisotropic convolution applied to a larger, circularly symmetric PSF. In this case, the 
effect of the tiny asymmetric deconvolution upon the weighted second moments of a given image 
can be expressed as a fourth-order weighted moment of the image, which KSB christen the smear 
polarizability. Given the smear polarizability of an image and a measure of the anisotropy of the 
PSF, the measured second moments are corrected analytically for the PSF bias. 

The FT97 method differs in that the correction for PSF anisotropy is applied to the image 
rather than to the measured moments: a minimal 3x3 convolution kernel is created which will 
"circularize" the PSF. The galaxy shapes are measured after this kernel is applied to the image. 

The primary drawback to these methods (K00) is that the approximation upon which they 
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are based often fails: a typical diffraction-limited PSF in no way resembles a small convolution to 
a round PSF, and even a simple aberration such as coma creates PSFs that violate this condition. 

These methods have features, however, that we wish to retain in any improved formulation. 
They are easily adapted to a PSF that varies across the image — the requisite moments of the 
PSF are measured wherever a star falls upon the image, and interpolated to the location of each 
galaxy. The FT97 method, by fixing the image, frees us from having to measure the higher-order 
moments that comprise the smear polarizability, though at the expense of a slight reduction in 
image resolution and/or an increase in the image noise. 

The perturbative methods correct only the PSF bias, not the dilution, because the dilution 
cannot be considered a small perturbation in any extant data set. FT97 calibrate the dilution with 
simulations, as mentioned above. The original KSB work made use of empirical calibration as well, 
as the shear polarizability (c/. §3.3) measures only the susceptibility of the image to a distortion 
which might be applied after the PSF convolution. Wilson et al. (1996) suggest an empirical 
calibration by deconvolving the real images, applying a shear, reconvolving, and remeasuring to 
determine the response. Luppino &; Kaiser (1997) introduce the pre-seeing shear polarizability, 
which approximates the susceptibility of the KSB weighted moments to a shear applied before the 
PSF. The pre-seeing shear polarizability is, to its lowest order, similar to the resolution factor R 
introduced above for unweighted moments. There are, however, fourth-order moments of galaxy 
and stellar shapes involved, and all are measured with Gaussian weights so the noise does not 
diverge. 

The KSB method updated to use the pre-seeing shear polarizability is exact for Gaussians and 
in the limit of a compact anisotropy kernel, but is not exact in the general case (K00). There are 
additional ambiguities regarding the appropriate size of the Gaussian weight to be applied when 
measuring the PSF moments (Hoekstra et al. 1998). The updated KSB method is in wide use, and 
several papers have investigated how accurately it performs for simulated galaxies and PSFs that 
are not Gaussian (Erben et al. 2001; Bacon et al. 2000). While it is clear that in many circumstances 
KSB is good enough, we would prefer to understand and overcome its limitations. 

We will demonstrate below that the KSB and FT97 bias corrections are the first terms in a 
series expansion of the deconvolution of the galaxy image. 

Rhodes, Refregier, & Groth (2000) investigate the KSB method in some detail, carrying forth 
the transformations equations to higher order than do KSB. The PSF corrections, however, still 
require substantial approximations. Below we construct a hierarchy of all the weighted moments of 
the image and give general formulae for their transformations under shear, convolution, and other 
operations. 
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6.1.4- Deconvolution Methods 

Kuijken (1999) suggests that one determine shear by summing the images of all the galaxies 
in some cell, then comparing this high-S/iV summed galaxy image to the PSF. The summed galaxy 
image should approach circular symmetry in the source plane; the shape of the summed image 
can thus be modeled relatively simply as a circular source with arbitrary radial profile, sheared by 
the lens to have elliptical isophotes, then convolved with the PSF. He adopts a flexible parametric 
representation of this mean radial profile. A candidate profile with candidate lens distortion can 
be convolved with the known PSF and compared to the measured mean image. The radial profile 
and distortion are then adjusted to give a best fit. This method can be viewed as a limited form of 
deconvolution: the only characteristic of the deconvolved image one cares about is the ellipticity; 
the multipoles beyond quadrupole are irrelevant to the measurement and are discarded. 

In Kuijken's method, the higher-order multipoles are discarded by averaging over many galax- 
ies. This may be difficult in practical situations, where the PSF and/or the distortion signal are 
varying across the field too rapidly to gather sufficient background galaxies to sum. Another draw- 
back is that summing galaxies' images may not be the optimal way to combine their information 
on the shear. Kuijken suggests applying the method to individual galaxies given these potential 
problems, which amounts to an assumption that the galaxy ellipticity is constant with radius. The 
accuracy of this method is not discussed in detail. But the method does have the strong advantage 
of being able to cope with arbitrary PSF behavior, and simultaneously removing both the bias and 
dilution effects from the measurement. 



6.1.5. Commutator Method 

K00 introduces a new approach to PSF correction, deriving from the PSF an operator which 
can be applied to the observed image to effect the transformation of applying a given pre-seeing 
shear. This operator is derived by considering the commutation of the shear and convolution 
operators in Fourier space. With the operator in hand, the response to a pre-seeing shear can be 
determined directly from the post-seeing image, given sufficiently detailed knowledge of the PSF. 

The K00 formulation is the first, to our knowledge, to offer an exact correction for PSF bias 
and dilution. We take a different approach below, and then offer some comparison of the two 
approaches. 



6.2. Optimal Methods in the Presence of a Convolution 

In §2.5 we noted that a real-space shear is equivalent to an opposing shear of the Fourier-space 
image, so that we can conduct the roundness test in Fourier space. One virtue of using a Gaussian 
weight function for our roundness test, as in Equation (3-21), is that this test takes the exact same 
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form in /c-space, namely 

M(I) = «=► = M(I) = !d 2 ke- k2a2 l 2 I{^)k 2 e- 2i ^. (6-3) 

The only difference is that the Gaussian weight has width \ja in /c-space, i.e. a broader weight 
in real space is narrower in Fourier space. Furthermore, our assumed uniform white noise in real 
space transforms to uniform white noise in Fourier space. Therefore the entire derivation of the 
optimal weight in §3.1 could have been executed in /c-space without any change in the result. 

The effect of the PSF convolution is to suppress the image by some transfer function T(k). 
With perfect knowledge of the PSF, we can deconvolve the observed image 1° to retrieve the 
image-plane transform P = I°/T (here we mean the image plane of the gravitational lens, before 
the application of seeing). A deconvolution would remove the PSF bias entirely, but the noise is no 
longer homogeneous in /c-space, having been amplified (perhaps infinitely) by 1/T. If the PSF is 
close to circularly symmetric, then T is nearly independent of direction. If we now imagine making 
our roundness test on the deconvolved /c-space image, we adapt Equations (3-7)-(3-9) to give 



SM = j / kdkk z w{k)k-± (6-4) 

Var(M) = —2 / kdkd<pk A w 2 (k) cos 2 2<p/T 2 (k). (6-5) 

^opt(fc) oc Z ^^- (6-6) 

The optimal filter is therefore narrower in /c-space than both dl/d(k 2 ) and T(k). Hence in real 
space, on the deconvolved image, the optimal filter is broader than both the object and the PSF. 
This means that we should, sensibly, restrict our roundness test to the region of /c-space that (a) 
has signal in the true galaxy image, and (b) is not suppressed below the noise by the convolution. 

Our strategy might then be to create some kind of deconvolved image, and apply a Gaussian 
weight to test for roundness in the deconvolved /c-space (which is equivalent to using a Gaussian 
weight in real space). For a Gaussian object with pre-seeing size it, and Gaussian seeing with 
size a*, the optimal weight in deconvolved k-space is a Gaussian with size (2a 2 + a 2 ) -1 ' 2 . Such a 
roundness test is equivalent to a roundness test on the observed, real-space image with a Gaussian 



weight of size \ af + a% = 0~ o . So the optimal size of the weight is again matched to the size of the 
observed image. 

Recall that our algorithm for measuring S requires that we shear the coordinates until the 
object appears round (M — ► 0). We want to apply this shear to the deconvolved image. If the 
object is not round to begin with, then in this sheared coordinate system the transfer function T(k) 
will no longer have azimuthal symmetry, which will invalidate the above derivation of the optimal 
weight. We will still have a valid measurement of the shape of the deconvolved object, but possibly 
with sub-optimal noise level. The increase in noise is a second-order effect, however, so we will not 
bother to re-optimize the roundness test for this asymmetric noise spectrum. 
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When the transfer function is anisotropic, then the noise spectrum in the deconvolved image is 
also anisotropic. There are subtle second-order effects, described in §8, which bias the orientation of 
the measured shear in the presence of an anisotropy in the noise. Such a noise anisotropy is present 
after any deconvolution of an anisotropic PSF. A noise anisotropy is also present in the observed 
image if the PSF is anisotropic and we are not strictly sky-limited. Furthermore, K00 points out 
a selection bias that can creep into the shear measurements even with a perfectly unbiased shape 
measurement algorithm. We will discuss means to defeat these biases in §8. 

Applying a Gaussian weight in deconvolved fc-space leads, formally, to divergent noise 'ifT(k) = 
for some k. This is a real problem, as any finite-sized telescope must produce a transfer function 
that is identically zero beyond some critical k c . The deconvolved image hence has infinite noise 
for k > k c , while the Gaussian weight remains finite. As we shear the fc-space to make our source 
appear round, the "wall" of infinite noise moves inwards to e~ r, ' 2 k c . In order for our method to 
remain feasible, our deconvolution algorithm must not attempt to fully deconvolve those portions 
of fc-space at or near k c , so that the noise remains bounded. There is hence a balance to be struck in 
executing the deconvolution: we want the carry out the deconvolution to sufficiently high order that 
the effects of the PSF upon the source ellipticity are removed, but we do not want to deconvolve 
high-order details that will increase the noise. It would seem, intuitively, that this is possible, since 
the ellipticity we seek, the Gaussian- weighted quadrupole moment, is a low-order characteristic of 
the deconvolved image. 

The method we describe below is based upon an expansion of the image and PSF into hier- 
archies of Gaussian-weighted moments — essentially an eigenfunction decomposition. Convolution 
corresponds to a matrix operation on the vector of moments. The moment vector, and hence the 
convolution matrix, are formally infinite, but we can choose to truncate the description at some 
order which we believe to contain all the useful information on the image ellipticity. The con- 
volution matrix is then finite and can be inverted, and the deconvolution executed as a matrix 
operation. The high-order moments are not deconvolved, so the noise in the deconvolved moments 
remains finite and in fact close to optimally small. This moment-based method reduces the de- 
convolution (as well as other transformations) to a matrix multiplication, which can be executed 
on an exposure-by-exposure basis even for objects with very low S/N on a single exposure, and is 
therefore very practical for the purpose of weak lensing measurements. 



6.3. The Laguerre Expansion 

6.3.1. Definitions 

The simplicity of the formulae for deconvolution in the special case of Gaussian objects, plus 
the utility of the Gaussian weight for shape measurements, lead us to seek a description of the 
image in some Gaussian-based expansion. To maintain the simplest form for convolutions, we look 
for a decomposition of our images into components which are eigenfunction of the Fourier operator. 
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Such functions will also be eigenfunctions of (—V 2 + r 2 ), hence we are led to the eigenfunctions of 
the 2d quantum harmonic oscillator (QHO). 2 In one dimension, the QHO eigenfunctions are each 
a Gaussian times a Hermite polynomial. The Edgeworth expansion, familiar to many astronomers, 
is a decomposition into Id QHO eigenfunctions. The 2d version we call the "Laguerre expansion," 
using the QHO eigenfunctions 

I(r,d) = J^M&M) ( 6 " 7 ) 

p,g>o 



r pq (r,0) = ^M^Te^e^l^L^^/^) (p > q) (6-8) 



{-l) q q\ (V 
\pKa 2 y p\ 

r qp ¥ m (6-9) 

m = p — q. (6-10) 

L q m '(x) are the Laguerre polynomials, which are defined by the generating function (Abramowitz 
& Stegun 1965) 

(1 -*)-«-! exp f -^- J = £ 2$ (*)*"*. (6-11) 

The Laguerre polynomials satisfy many recurrence relations; the following provide a way to calculate 
them rapidly: 

4 m) (x) = 1 (6-12) 

LJ m) (x) = (m+l)-x (6-13) 

(q+l)L { ™\(x) = [{2q + m + l)-x\L q m \x)-{q + m)L q ^\(x). (6-14) 

A QHO with wavefunction ifj pq has angular momentum m = p — q. We will also make use of the 
quantum number N = p + q, which is the excitation energy of the state. Any two of {N, m,p, q} 
suffice to specify the state. The intensity multipole functions I m (r) are composed from the ipNm 
for N = \m\, \m\ + 2, \m\ + 4, . . .. The polynomial in ip pq has terms up to order N in r, and ip pq 
can also be expressed as the Gaussian times a (complex) polynomial of order N in x and y. A few 
low-order ^ vq are plotted in Figure 6. 

Refregier (2001) has independently introduced the application of QHO eigenfunctions to galaxy 
shape analysis. Some of the results and ideas presented in this section are presented therein, and 
applied in Refregier & Bacon (2001), with different notation. Refregier (2001) also presents useful 
relations for a Cartesian-based family of 2d QHO eigenfunctions. We make use only of the polar 
family, which are eigenstates of the angular momentum and hence have strong rotational symmetry. 

The Laguerre functions have many properties we will find useful. As eigenfunctions of the 
harmonic oscillator, they are orthonormal, up to the factor a which we introduce in order to give 



Thanks again to the anonymous referee for providing this logic. 
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Fig. 6. — The first few of the orthogonal functions ip pq are plotted. Only the real parts are plotted, with 
blue (red) for positive (negative) regions. The characteristics are familiar from their use as eigenfunctions 
of the 2d quantum harmonic oscillator: each has 2m = 2\p — q\ azimuthal nodes and N — 2m radial modes, 
where N — p + q is the number of quanta in the state. The i/'20 component is most important as it responds 
in first order to shear, and hence its absence is our test for "roundness." 
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the b pq units of flux: 



A -. - 1 



dxtpp q (y:)ip^ q ,(x.) = —5 pp >S qq > (6-15) 



a 



b pq = a 2 /d 2 x/(x)^(x). (6-16) 



Thus b pq is a Gaussian-weighted moment of the intensity. Since / is real, we must have b pq = 
b qp - When measuring the b pq from an image with uniform white noise, Equation (6-16) yields a 
covariance matrix for the b pq that is diagonal: 

Cov(bpqb p > q/ ) =na 2 5 pp >5qq>, (6-17) 

where n is the number of counts per unit area. The variance in b pq is shared between the real and 
imaginary parts, except for b pp , which must be real. Non-uniform noise, e.g. shot noise from the 
galaxy itself, produces a more complicated covariance matrix, as described in Refregier (2001). 
The significance v of detection with the Gaussian filter is given by 

v 2 = ■%. (6-18) 

Our algorithm for measuring object shapes requires finding the coordinate system in which 
the centroid is zero (Equation [3-18]) and the roundness test (Equation [3-6]) yields zero. When 
disregarding seeing, we also set the weight size a by maximizing the significance (Equation [3-19]). 
These conditions are succinctly stated by the Laguerre coefficients: 

Centroid: b 10 = (6-19) 

Roundness: 620 = (6-20) 

Significance: 611 = 0. (6-21) 

The first two of these equations involve both real and complex components, the third is real. We 
must satisfy these equations by translation, shear, and dilation of the object (or of the underlying 
coordinate system). These operations can be expressed as transformation matrices acting upon 
the vector b = {b pq }. The determination of shape is thus efficiently executed by measuring b in 
the original coordinate frame, converting Equation (6-16) to a sum over pixels. Henceforth we can 
iterate to a solution of our three conditions by manipulating b, and there is no need to return to 
the pixel data. 

In the case of finite resolution, we wish to satisfy Equations (6-19) and (6-20) for the decon- 
volved image. We can express the deconvolution as a matrix operation on b as well. So we need 
to find the matrix equivalents of the translation, dilation, shear, and convolution transformations. 
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6.3.2. Raising and Lowering Operators, Point Transformations 

The raising and lowering operators for the 2d harmonic oscillator wavefunctions have the 
properties: 



a; = 
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a p ' = 
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a° = 
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17 I dx Z 9?; 



"■p r pq 
a P Vpq 

a q 1 Ppq 
a q Wpq 



VP^(p-l)qi 
\ / P+ T ^(p+l)q^ 

v^V>p(,_i), 

v^+T^pVi)- 



(6-22) 



Note that the commutators of the operators are all zero except for [a p ,a P ] = [a,j,aj] = 1. These 
operators can be used to formulate the transformation matrices we need. Consider first a dilation 
of the image / by factor 1 + \x with ^ <C 1 to a new image V: 



I'(x,y) 



J[(l - n)x, (1 - n)y] 



, d d 

ox oy 



H 



pq 



(6-23) 
(6-24) 

{l-p(aX-«-l)}j (6-25) 

^*w{(l + /j)i|' M -/iv , M^-i)(,-i) + IJ-y/ip + l)(g + 1 ip( P +i)(q+i) } (6-26) 
(1 + n)b pq - /V(p + l)(q + 1) 6(p+i)( g+ i) + /ivS^p-iK?-!) (6-27) 



An infinitesimal dilation |U thus mixes coefficients (N ± 2, m) into the (N, m) coefficient — as ex- 
pected, multipole order m is conserved by the transformation. The dilation transformations arc 
thus generated by the matrix 

d = 1 - a u v a u q + apa^, (6-28) 

and the transformation matrix D for a finite dilation e^ can be expressed as 



D, 



exp [fid] . 



(6-29) 



The finite dilation will preserve m as the generator does, but will mix terms of N + 2j together, 
with j any integer. The matrix elements for the finite dilation are derivable in closed form and in 
a rapid recursion form; the latter are derived in Appendix A. 

The generator s for the shear transformation is also easily expressed in terms of the raising 
and lowering operators. For a shear ?j<1 along the x axis, 



I'(x,y) - I[(l- V /2)x,(l+ri/2)y] 
-li^x-^y)} 1 



;i + 7?s)/, 



i 



(aj,) 2 + {a\f -a 2 p -a 



(6-30) 
(6-31) 
(6-32) 
(6-33) 
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=* K, = ^ + l(Vp(f^e~ 2l %-2) q + Vq(^e 2lf3 b p(q ^ ) 

-V(P+l)CP + 2)e 2 % +% - V(9 + !)(? + 2) e- 2i % (9+2) ) , (6-34) 

S^ = expfrys]. (6-35) 

In (6-34) we have inserted the phase factors which result from a shear at arbitrary position angle 
[3. To leading order, the shear mixes in states with p or q changed by ±2; in the finite case, the 
shear transformation mixes together b pq for which p and q each change by any even number. Note 
that this equation generalizes the weighted-moment transformations in Rhodes, Refregier, & Groth 
(2000) to all orders. Note also that the KSB method identifies the quantity 620/^00 as the ellipticity, 
so the "shear polarizability" is compactly expressed using Equation (6-34) with p = 2, q = 0. 

An infinitesimal translation of the image by (xo,yo) gives 

l'(x,y) = I(x-x ,y -yo) (6-36) 

1 -^i--yoi-)i (6-37) 



dx dy / 

1 + zt- ztf)I, (6-38) 



xo + iyo 



(6-39) 



z 



=> b 'pq = b PQ + g {Vq~ b p(q-V ~ VP + 1 b ( P +l)q) + g (v^ 6 (p-l)g ~ V Q + 1 fe p(<?+l) j ( 6 " 41 ) 

T 2 = exp[zt-zt f ] (6-42) 

To leading order the translation changes p or q by ±1, and the finite translation T z mixes all the 
states together. Appendix A gives the matrix elements for all the finite transformations. 

The "smear polarizability" coefficients of KSB could be quickly derived at this point by noting 
that an infinitesimal convolution along the x direction can be expressed as an average of two images 
translated by ±xo, which will lead to a smear generator that is second order in the raising and 
lowering operators. Appendix C gives a closely related derivation. 



6.3.3. Shape Measurement from Laguerre Decomposition 

Our algorithm for shape measurements (with perfect resolution) is to find the translation z, 
dilation fj,, and shear rj which must be applied to the image to satisfy the conditions (6-19)-(6-21). 
If b is the vector of Laguerre coefficients for the image I, we can write 

b' = M ■ b = (S^D^T,) • b, (6-43) 
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and the elements of b' with (pq) = (10), (11), and (20) must vanish. This is in general a complex 
non-linear equation, but for /i, 77, z <C 1 the linearized Equations (6-27), (6-34), and (6-41) yield 

&io = &io + (&io- V2b 21 )^-^b 12 V-^hoV 

-^hoz + ^boo-b^z (6-44) 

b' u = b 11 + (boo + b 11 -2b22)l^-^b 13 fj-^b 31 ri 

+i(6io - V26 2 i)^ + 5(601 - v 7 ^)^ (6-45) 

6' 20 = 6 20 + (620 -V363i)m + ^(600 -b 22 )fj-&b m n 

-^b 30 z + i(v^6io - 621)2, (6-46) 

writing 77 = 77+ + ir/ x • These linear equations may be solved explicitly for the desired transformation 
coefficients 77, /i, and z. The solution appears much simpler if the object is nearly round, such that 
b pq <C 600 f° r P ¥" Q- I n this case the transformation parameters are, to leading order, 

- - 2 ^» (6-47) 



600 — 622 
-2601 



600 -6 



11 



(6-48) 



" = i bl L ■ ( 6 - 49 ) 

000 — -^22 

From these equations it is clear that 620, 610, and &n are the primary carriers of information on 
shape, centroid, and size, respectively. Since the shape component 77+ is (77 + fj)/2, its uncertainty 
is to leading order 



a = Var(r? + ) = — ^ (6-50) 

' (600 - b 22 ) 2 

Ana 2 

= ~w 



2 _ , r , , _ 2 [Var(6 20 ) + Cov(b 2 o6o2) + Var(6 02 )] 

-h 

[1 - 6 22 /6oor 2 (6-51) 

=> v v = - [1 - W&ooP 1 • (6-52) 

v 

Here we have made use of Equation (6-17) for the covariance matrix in the case of white noise, 
and Equation (6-18) for the definition of the significance v. We see that the Laguerre expansion 
easily reproduces the earlier result in Equation (3-26), with the 04 parameter in that equation being 
simply the strength of the b 22 coefficient. 

With a more tedious procedure we may derive the solution for 77 to second order in b pq /boo. 
We take b pp <C 600 for p > this time in order to simplify the expression: 

^ - 2 ^^fe-2v / 3^ 

, 2tf n -2\/6biabo3-2\/2b ibi2 (6-5^ 
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The terms in the first line arise from solving for 77 = 0; the terms in the second and third lines 
arise when simultaneously solving the centroid and size constraints, respectively. From these terms 
one can derive the 0(v~ A ) terms in the shape uncertainty if desired. But more importantly, this 
expression uncovers a very significant source of bias in shape measurements: the presence of the b pq 
in second order means that it is possible for noise to be rectified in the determination of e. If the 
noise is anisotropic — as is the case after correction for an anisotropic PSF — then e will be biased. 
We will explore this in more detail in §8. 



6.3.4- Fourier Transforms and Convolution 

The observed galaxy intensity /°(x) is the convolution of the image-plane intensity -P(x) and 
the stellar PSF J*(x): 

J°(x) = r(x) o I*(x) = jd 2 x' /* (x')/*(x - x). (6-54) 

The post-seeing, pre-seeing, and stellar images can be expressed as vectors b°, b 4 , and b*, respec- 
tively, of coefficients over our eigenfunction sets ipZs, ippL and ipZq. The convolution can then be 
expressed as a matrix relation 

b° = C(b*)-b* (6-55) 

b U = E C plT q %m b U- ( 6 - 56 ) 

We will effect the deconvolution by inverting a truncated version of the matrix C. Its coefficients 
are determined by the relation 

r P ; qi o ru = E cv vX^%u - ( 6 - 57 ) 

The convolution is more easily expressed in /c-space. With the usual definitions (denoted as "System 
3" by Bracewell (1978)) 



J(k) = (27T)- 1 d 2 xl{x)e- ik - x , (6-58) 

I(x) = (27T)- 1 fd 2 kl(k)e ik - x , 

the convolution becomes a multiplication, and the matrix coefficients in Equation (6-57) can also 
be expressed as 

2^fefe* = E C fX^lW ( 6 - 59 ) 

We now make use of another remarkable property of the Laguerre-exponential eigenfunctions, which 
is that they are their own Fourier transforms. First we note that ipQ is a two-dimensional Gaussian, 
and the transform (6-58) is easily executed to yield a new Gaussian, 

Co = -L e - fc ^/2. (6.60) 

\/7T 
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The functional form of ipoo(k) matches that of ipoo(x), except that we must send a — > 1/cr. Next 
we can use the definition of the raising operators (Equation [6-22]) to note that 



°1 / 

a p ] ip 



(7 ( Kx ~\~ 1 f^xi 



l 



a \dk x 



d d 

+ i- 



dh. 



4> = ~a a p H 



(6-61) 



Thus the fc-space raising operator has the same form as the x-space raising operator, save the 
a — ► 1/cr transformation and an additional factor of —i. The same is true of a q . Since ipQ$ and the 
raising operators are each unchanged by the Fourier transform, it must be true that all the ip pq are 
their own Fourier transforms as well, with appropriate scaling factors of i and a. More precisely, 
we have (c/. [6-8]) 



#,(M) = { S- M {ka) m e m <t>e- k2 ° 2 l 2 L q m Xk 2 o 2 ) (p > q) 



(6-62) 



So the problem of convolving two of our eigenfunctions is reduced to the simpler problem of mul- 
tiplying the same two eigenfunctions. We may reach several conclusions immediately: 

• The convolution of V'^. m . with tpff will produce an observed image with azimuthal order 
m = rrii + m*. Furthermore the multipole phase of the observed image must be the sum of 
the phases for the original and PSF images' coefficients. 



If we choose the scale size a for the observed image eigenfunctions such that 



o-: 



af + at 



(6-63) 



then the convolution contains components ip^f only for N < Ni + N±. Recall also that we 
must have iV > m = rrii + m*, so it must be true that ipp* o ipplq ^ ^v+v )a ^ or 1 = ®- 

In Appendix B we give a recursion relation to calculate any of the elements Cplqf* q * that we 
need to calculate the convolution, and thus the deconvolution. 



6.3.5. Deconvolution, Noise Amplification, and Truncation 



With the formulae for the convolution matrix C in hand, we can investigate the nature of the 
trade-off between the fidelity of the deconvolution — i.e. the extent to which effects of the PSF upon 
the shape are removed — and the noise level of the deconvolution. We must truncate the b vectors at 
some finite order N to implement the deconvolution; higher N will remove more of the PSF effects, 
but also increase the noise level. We can illustrate this phenomenon by considering the simplest 
case of convolution by a unit Gaussian PSF, 6* = 2^/^T5pQ5 q Q. Since the PSF has only m* = 
terms, we will have m = ?rtj, so the convolution matrix is block-diagonal and we may deconvolve 
each m independently. If we choose o~1 = o~f + a+, then the convolution matrix coefficients are zero 
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for N > Ni, and hence each block matrix is also upper-triangular. This makes the inversion easy. 

Using the results of Appendix B, defining a deconvolution parameter D = crf/a^ = 1 — crl/a^, and 

using Equation (6-17) for the covariance of the measured moments, the deconvolved value of 6q 

and its variance are 

00 /i _ n\P 

6 oo = E(-d-) (-!)%> t 6 " 64 ) 

p=0 ^ ' 

=* Var(6* 00 ) = n,^ -^ = n *%TT (^ > V2)- (6-65) 

In practice we could not measure the b° to infinite p and we must truncate the sum (6-64) at some 
finite p. The more terms we include, the more accurately we describe the deconvolved 6q , but each 
additional term in the deconvolution adds more noise. For D > 1/2, the added noise in successive 
terms drops as p increases, and we could in principle do a complete deconvolution of b l 0Q with finite 
variance. For D < 1/2, however, the noise in the deconvolved moment diverges for p — ► oo, and we 
are forced to truncate the deconvolution matrix. 

For ellipticity measurements we are primarily interested in the 620 moment. For this case of a 
Gaussian PSF, the general b^ moment deconvolution is 

00 

";..o = Et^i^i i-D'i/r:"]^), ^ 




***-> - ^E(^rcr) <«*> 



<?=0 



j~, \ m+] 



= na i D [2D^l) ( jD>1 /2). (6-68) 

We note that these moments also demonstrate the properties that they are noisier than the observed 
moments, and infinitely so if D < 1/2 and the deconvolution is not truncated. The noise increases 
as m increases, as we would expect, since higher spatial frequencies must be recovered. 

For a purely Gaussian PSF with D > 1/2, it is possible to complete the deconvolution with 
finite noise. A real PSF must have components beyond 6g , since the transfer function vanishes 
above a critical k value. When we invert the convolution matrix for such a PSF we will find that 
the variance series akin to (6-65) would diverge for any D < 1. The best truncation value for the 
matrix will depend upon the form of the PSF and the accuracy to which we demand correction of 
PSF effects. We will examine some specific cases in Jarvis et al. (2002). 



6.3.6. Analytic or Kernel Correction for PSF Bias? 

There are two general means of eliminating the shape bias induced by the PSF. One alternative 
is to measure the anisotropy of the PSF carefully and apply analytic corrections to the measured 
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objects, as occurs naturally within the deconvolution framework described in the preceding para- 
graphs. The KSB formalism contains an approximate analytic correction. The other method is to 
convolve the image with a spatially-varying kernel which removes the anisotropy from the PSF, as 
first demonstrated by Fischer & Tyson (1997) and further advocated in K00. Removal of the PSF 
bias is the most critical task for a weak lensing method: the PSF dilution is only a calibration error 
or signal modulation, but the PSF bias introduces a first-order false signal to the lensing analysis. 

Ideally, analytic correction is preferred. The convolution kernel will always degrade the im- 
age resolution to some extent. Many PSFs cannot be "rounded out" unless the kernel is of size 
comparable to the PSF itself (K00). Furthermore the convolution will perturb the noise spectrum 
of the image, complicating error estimation. The kernel method, however, can for simple kernels 
be faster than an involved analytic correction, especially if eliminating the bias is more important 
than calibrating the dilution for the task at hand. 

The kernel method may offer peace of mind. The demands for rejection of the PSF bias are 
very stringent: PSF ellipticities of 10% are not uncommon, yet state-of-the-art lensing surveys 
require systematic errors below 0.1%, so the method we choose must reduce the PSF bias by 
over 2 orders of magnitude. If the kernel method can make the PSF truly round, then symmetry 
principles preclude any artificial coherent ellipticity. We flesh out this statement in §7. We may 
have more trust in the symmetry principle than we do in an analytic correction, especially when 
finite sampling is taken into account. Note, however, that the symmetry principle requires that the 
noise power spectrum is also isotropic; §8 below demonstrates how anisotropic noise can bias shape 
measurements through centroid bias. 

We therefore will describe an alternative procedure: eliminate the bias with a kernel correc- 
tion; measure the shapes, either isotropizing the noise or making a correction for centroid bias; 
then correct them for dilution with a heuristic formula such as Equation (6-1). In §7 we describe 
a rounding-kernel methodology, and in Appendix C we derive a higher-order version of the dilu- 
tion correction. These methods were used by Fischer et al. (2000), Smith et al. (2000), and the 
forthcoming CTIO lens survey. 



6.4. Pixellated Data 

Real astronomical data is integrated into pixels and sampled at finite intervals, and the con- 
tinuum limit that we have assumed in our analysis is not strictly valid. A formal description of 
the process is that the PSF is altered by convolution with the pixel response function (PRF) to 
produce the effective PSF (ePSF), which is then sampled at the pixel pitch a — see Lauer (1999) or 
Bernstein (2001) for further elaboration. If the pixel pitch is coarser than the Nyquist interval of 
the PSF (X/2D for a diffraction- limited image), then there is aliasing and we cannot unambiguously 
reconstruct the original image or the true b pq coefficients. 

In the case where the pixel size a is small compared to the PSF size, the b pq may be estimated by 
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converting the integral Equation (6-16) into a sum over pixels. But the formally correct procedure 
is to fit the vector of sampled pixel fluxes I{ with a model galaxy with free b° . This is a standard 
minimization problem for the \ 2 parameter 






2 ^-E^AV^ 

v ^ — ^m — (6-69) 



Since the model is linear, the solution for b° and its covariance matrix have a closed form, but of 
course the number of b° coefficients allowed in the model must be less than the number of pixels 
being fit. With the best-fit coefficients and their covariance matrix in hand, we may proceed with 
the methods outlined above. It is the ePSF rather than the PSF that has convolved the source 
image, but since we in fact measure the ePSF for stars, the deconvolution procedures are unchanged. 
Dithered exposures can be handled by extending the pixel sum over multiple exposures, as long as 
the PSF is unchanging. Poorly sampled images will manifest the aliasing as strong degeneracies in 
the solution for the b pq . 

A potential time-saving procedure would be to fit the pixel data to a Laguerre expansion that 
is already convolved with the local PSF, minimizing 
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(6-70) 



Here C is the Laguerre-coefficient convolution matrix for the local PSF. The intrinsic coefficients 
b % are the unknowns in this linear fit. This is very reminiscent of the method of Kuijken (1999), 
except that our model for the intrinsic galaxy shape may depart from circular symmetry, and our 
Laguerre formalism will allow the fit to proceed as an matrix solution, with matrices rapidly built 
by recursion. Aside from its suitably to poorly-sampled data, this direct-fitting approach has two 
further advantages over the deconvolution-matrix method of the previous section. First, when 
the noise is not flat and Equation (6-17) does not apply, the direct-fitting method more directly 
produces the full covariance matrix for the b* vector. Second, the direct fit may make use of images 
that are partially contaminated by invalid pixels, e.g. cosmic rays, since they may be excluded 
from the \ 2 sum. It is this method which we consider most promising for real data. 



7. Rounding-Kernel Methods 

In this section we develop a method for producing a convolution kernel that symmetrizes the 
PSF to some desired degree. As discussed above, this is a potentially efficient means of reduc- 
ing or eliminating PSF bias from the shear determination. By taking advantage of the Laguerre 
decomposition, the derivation and application of the spatially-dependent kernel can be efficiently 
implemented. 
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7.1. The Target Transfer Function 

As discussed in §6.2, the effect of seeing is to convolve the initial image, P, with a PSF transfer 
function T to produce an observed image, 1°. One could theoretically remove all the effects of seeing 
by convolving the observed image with a kernel K whose Fourier transform K = 1/T to produce a 
final image, I' . 

However, as discussed above, this is impossible to carry out in practice, because T(k) = for 
all k above some critical value. We can avoid this problem by ignoring the high order moments of 1° 
and P . This is straightforward when one does an eigenfunction expansion of the images using the 
eigenfunctions introduced in §6.3, since the expansion can be truncated at some value of N = p + q. 
This captures most of the real information about the PSF without introducing the noise from the 
higher order moments. After the convolution, we then have a target transfer function, T", which is 
not exactly a delta function, but which can be made to match a delta function up to some order 
N. Most generally we have 

1° = Top (7-1) 

pq 

if = T' op = KoToP (7-3) 

T ' = £OwM) ( 7 - 4 ) 

pq 



For the ideal case of T' = <5(x) 



(-l)P 



Thus, if X" satisfies Equation (7-5) up to some cutoff order, TV, I* will be approximately identical 
to P up to the same order. If we can find a kernel K, which produces this target transfer function, 
T', we will be able to remove all the effects of seeing as well as possible given our ignorance of the 
high order terms. 

How stringently do we need to satisfy Equation (7-5)? Less strict requirements on T" make 
it easier to find an appropriate, compact kernel. Our present goal is to create a transfer function 
T' which does not produce any shear bias. If the original scene /' is unlensed, and we represent 
the shear measurement process as some operator 5(1), then the isotropy of the Universe guarantees 
that S(P) = (we consider S as a complex number 8+ + i8 x )- Our demand on the transfer function 
is that 

5(lf) = 5{T' o P) = 0. (7-6) 

Most generally the results of the shear measurement process can be expanded as a power series in 
the coefficients of the T' PSF: 

0(1 O 1 ) = 2_^ a pi q 1 piqi + 2_^ a Piqi,P2q2"piqi< ) p2q2 + " " " (.'"') 

piqi pi<2i:P2<?2 
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where the en are some coefficients that depend upon the measurement process, the image charac- 
teristics and the size of the PSF. We now examine the consequence of rotating the image and the 
PSF by some angle (3. The measured shear and the T" coefficients behave as 

5(I f ) -► e 2 ^5(I f ); &£ -► e { ^^b% q . (7-8) 

The individual galaxies in the original scene I* are not invariant under rotation, but any statistical 
measure of their collective properties must be invariant under rotation, i.e. P is invariant under 
rotation in the same sense that the Universe is isotropic. Therefore the a% coefficients of Equation (7- 
7) are unaffected by the rotation, and to satisfy the conditions (7-8) we must have 

a piqi = for pi-qi^ 2; 

for pi - q 1 + p 2 - g 2 ¥= 2; (7.9) 



Thus only PSF terms with m = p — q = 2 can cause a shear bias, to first order. The primary goal 
of our kernel, therefore, will be to set 

b* p ' q = (m = p - q = 2) (7-10) 

With this condition satisfied, shear bias can only be of order (b*) 2 where p 7^ q. Higher-order a's 
are non-zero only if ra\ + m<i + • • • = 2. Many of these elements are zero as well; for example by 
considering the invariance of 5(1) under infinitesimal translation, one can demonstrate that aio,2i 
vanishes, but aio,io must exist. Satisfying (7-10) does not, therefore, guarantee the elimination of 
shear bias to all orders. 

For absolute assurance that a shear bias is absent, we can set 

b;' q = (p - q + [mod 4]), (7-11) 

i.e. insure that the PSF has fourfold symmetry. In this case the condition Y2 rm = 2 can never 
be satisfied for non- vanishing coefficients of T". Note that it would also suffice to enforce any m- 
fold rotational symmetry on the PSF beyond m = 2, e.g. the diffraction pattern of a triangular 
secondary-support structure cannot by itself cause shear bias. 

In practice we produce a kernel K which enforces condition (7-10) up to some order N . We 
can set b^ = by appropriate choice of PSF center. The remaining shear biases must be of order 
(b2i/boo) 2 , ^4i&oi/^oo> e ^ c -' which are generally quite small. 



7.2. The Components of the Kernel 

Since we are going to convolve the kernel with a pixellated image, we must construct the 
kernel as a 2d array instead of a continuous function. The simplest kernel is the identity kernel, 
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composed of an array of O's with a single 1 in the middle. This kernel conserves flux; we want a 
flux-conserving kernels, so we only consider kernels which are the identity kernel plus a kernel for 
which the elements sum to zero. 

Next note that taking the derivative of an image is usually approximated by a discrete differ- 
ence. For example: 

dl (I(x + dx) — I(x — dx)) 

dx 2dx 

Another way to write this equation is as a convolution: 



(7-12) 



Similarly, 



Ox 



m 

dy 





-1/2 1/2 | ol (7-13) 





(7-14) 




In fact, all partial derivatives of any order in x and y can be represented as a convolution 
by a zero-sum kernel. First and second order derivatives can be effected using 3x3 kernels; third 
and fourth order derivatives require 5x5 kernels, and so on. One can therefore think of the kernel 
as being made up of a sum of these derivatives (d / dx) 1 (d / dy) 3 including, of course, the identity 
kernel, (d / dx)° (d / dy)° . 

The eigenfunction expansion of our images actually suggests a slightly different set of compo- 
nents for the kernel. Namely, 

K = Y J h j D ij (7-15) 



n - (— —X(—-—Y 
13 \dx dy J \dx dy J 

= jn{<-<y{ a p-<y ( 7 - 17 ) 

These components make it easy to use the raising and lowering operators to determine how a 
given kernel will act on an image. Note, however, that Dij are complex, which means that to end 
up with a real image after the convolution, we require that kji = ky. 

The 3x3 kernel components are: 

Dio = I -1/2 1/2 1 (7-1.8) 
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(7-19) 
(7-20) 
(7-21) 
(7-22) 



We note at this point that the FT97 method is equivalent to use of just the D20 an d -D02 kernel 
elements. The 5x5 and larger kernel components can be found by convolving the 3x3 components. 
For instance, D30 = D±q o D20, and D41 = D\\ o D 30 . 



7.3. Calculating the Kernel 

7.3.1. The Kernel for Infinitesimal Pixels 

The kernels given by Equation (7-18) - Equation (7-22) are only equal to the continuous 
derivatives (Equation (7-16)) to lowest order in 1/cr*. Typically, the PSF is only a few pixels in 
size, so this is not that good an approximation. However, it is a good place to start, as most of the 
technique will apply to the case of finite pixels. 

Combining Equation (7-15) with the definition of T" (Equation [7-4]), 

b*' = Y^ kijVijb* (7-23) 

ij 

The operator matrix Dy, in the limit of infinitesimal pixels, is easily calculated using Equa- 
tions (7-17) and (6-22), and there is a fast recursion 

D 00 b* = b* (7-24) 

D (i+1)j -b* = ^(<-<t )(D . ib *) 

Since b* is measured and the D^ are fixed matrices, we have a matrix equation for the unknown 
kernel coefficients kij, given the chosen constraints on b*' (e.g. Equation [7-10]): 



Mk = b +/ (7-25) 

where k = {kij} and the ij row of M is given by (Djjb*) T . Thus, given b* and b*', one calculates 
M using Equations (7-24), and then simply solves Equation (7-25) for k. 
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The problem with this method is that b*' will not usually be completely specified. Neither 
Equation (7-10) nor (7-11) fully constrains b*'. The easiest way to deal with this is to set k pq = 
for each unspecified b*' . Then the kernel will be as simple as possible while still satisfying all of 
the requirements for b*'. A more sophisticated technique for dealing with this issue is described 
below in §7.4. 



1.3.2. The Kernel with Pixelization 

The same general technique applies to the case of finite pixels in that we want to solve Equa- 
tion (7-25) for k. The only difference is the calculation of Djjb*. In this case, we must use 
Equations (7-18) - (7-22) rather than Equation (7-17). 

Consider the version of D\q given in Equation (7-18): 

D 10 b*(x,y) = i(b*(x+l,j/)-b*(x-l,j/)) 

+l(b*(x,i/ + l)-b*(x,i/-l)) (7-26) 

= i (T 2l b* - T_ 2l b*) + % - (T z2 b* - T_ 22 b*) (7-27) 

where zl = 1/c, zl = iju and T z is defined in §A.l. 
Calculations for the other Dy are similar. 



7. 3. 3. The Kernel in a Distorted Frame of Reference 

Most telescopes, especially those with large fields of view, have fairly significant distortion. 
To deal with this correctly, the shape measurements should be made in the undistorted world 
coordinates rather than the chip coordinates. The kernel's pixel grid, however, is still is in the 
distorted coordinate system. The calculation of Djjb* must therefore take into account the two 
different coordinate systems. 

If the world coordinates are (u, v) and the pixel coordinates are (x, y) then the correct values 
for zl and z2 are: 



\dx dx ; 

1 f du dv\ 

z2 = Aei + %) 
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7.4. Minimizing the Noise from the Kernel Convolution 

The motivation for having a compact kernel is to minimize the noise added by the convolution. 
A large kernel will use data with significant noise but little signal, adding to the noise in the 
convolved image. 

Therefore, let's consider the noise in the convolved image. Define K(m, n) to be the total 
convolution mask. 

K(m,n) = y^kijDij(m,n) (7-30) 



/•' (x, y) = 2_, K(m, n)I°(x + m,y + n) (7-31) 

mn 

If the noise in the image is dominated by sky noise, then we can define n 2 to be the variance 
in each pixel of the observed image. The noise in each pixel of the convolved image is 



n) = Y^ K(m, nfn 2 = n 2 ^ [ ^ k^D^m, n) ] (7-32) 

In §7.3.1, we set some kernel coefficients to to account for the unspecified components of b*'. 
Really, each of these coefficients is arbitrary, so the solution to Equation (7-25) will be 

k = k + AY (7-33) 

where ko is a specific solution, A is a matrix, and Y is an arbitrary vector. The dimension of Y 
is equal to the number of unspecified components of b*'. 

We can then find the particular Y which minimizes the noise added to the image by the 
convolution. For this derivation D(m, n) is the vector of each Dy (m, n). So, K(m, n) = D(m, n)k. 



dn 2 t 

(7-34) 



dY 

1 ^(D(m,n)kr 



d 



dY 

\ mn / 

= 2n 2 J2(D(m,n)k)(D(m,n)A) 
mn 

= ^(D(m,n)(k + AY)(D(m,n)A) (7-35) 
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Y^ (D(m, n)A) (D(m, n)AY) = - ^ (D(m, n)A) (D(m, n)k ) (7-36) 

mn mn 

This is now a matrix equation which can be solved for Y, which then gives you the solution 
for k which minimizes the noise. 

It turns out that if this method is implemented exactly as described, one ends up with a fairly 
large kernel, which is essentially a smoothing filter. While this solution will work, it is not exactly 
what we want from the kernel. We would rather have a smaller kernel which minimally changes 
the image. So, rather than minimize the noise in J* , we minimize the noise in I* — 1°. In other 
words, we leave out the Doo term in the vector products of Equation (7-36). This results in fairly 
compact kernels which vary smoothly across the image. 



7.4-1- Additional Kernel Components 

There are 9 free parameters in a general real 3x3 kernel. Equations (7-18) - (7-22) define 5 
kernels. The identity kernel is another. Thus, there are 3 more independent 3x3 kernels we could 
construct. Two of these can be made to approximate discrete versions of the derivatives D\q and 
Z?oi- The other can be made to approximate a discrete version of D\\. These alternate versions 
are: 



■l + t)/4 (l + i)/4 
AltDiQ - \ | (7-37) 

(— 1 — z)/4 (1-0/4 

(7-38) 
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(7-39) 



We can construct alternate higher order kernels in the same way as we constructed the regular 
higher order kernels. Namely, AUD^q = AltDio o D20, AUD41 = AltDn o D30, etc. 

These extra kernel components are useful if one is minimizing the noise as described above, 
since they add extra degrees of freedom for the minimization, and therefore can result in a smaller, 
less noisy kernel. 



7.5. Interpolating Across an Image 

The above discussion explains how to find the appropriate kernel given a particular PSF. 
However, in a real image, the PSF varies across the chip. Thus, the kernel will also vary across the 
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chip. 

There are two potential ways to deal with this. One can find the appropriate kernel for each 
star in the image. Then fit the kernel components kij as functions of (x,y). Alternatively, one can 
fit the coefficients b* of the measured PSF decomposition as functions of (x, y) and then solve for 
the appropriate kernel at each point. 

We choose the first method in our analysis for two reasons. First, the derived kernel can be 
directly applied to each star to make sure that it really does make the star round. Occasionally, 
a star will have significant high order components due to crowding or an uncorrected cosmic ray. 
When this happens, the derivation above fails, and we reject this kernel from the fit. It is a little 
cleaner to recognize these outliers with kernel interpolation than with PSF interpolation. 

The second reason to prefer fitting the kernel rather than the PSF is computational efficiency. 
The appropriate kernel must be calculated at pixel. It is significantly faster to evaluate a function 
than to solve a matrix equation. On a 2K x 4K chip, there are 8 million kernel evaluations. Both 
methods gain significantly by using a locally linear approximation to the spatial variation, but there 
is still a significant difference in computation for the two methods. 

The kernel scheme described here allows an substantial speed gain over a typical convolution 
method. Fourier methods are fastest for large convolutions, but are not practical for spatially- 
varying kernels. In our scheme, the output image can be written as 

I f = Y t k ij (x,y)*(D ij I°). (7-40) 

ij 

Instead of calculating the kernel at each output pixel, we can instead produce the images Dijl° 
(which are constant-kernel convolutions), and then accumulate sums of these images with spatially 
varying weights kij(x,y). More importantly, the Dijl° can be produced by recursive application of 
the 3x3 kernels. Hence the convolution by a 7 x 7 kernel can, for example, be reduced to three 
successive applications of 3 x 3 kernels. 

Note that interpolation of the PSF elements across the image is required for the analytic 
deconvolution described in §6.3.5. 



7.6. Dilution Correction 

Once a kernel has been applied to the image to symmetrize the PSF, the measured galaxy 
cllipticities need to be scaled by some resolution factor R to account for the PSF dilution. Ap- 
pendix C describes the scheme we use for estimating R, which closely resembles the KSB "smear 
polarizability" derivation. 
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8. Centroid and Selection Biases 

Even if the asymmetries of the PSF have been perfectly removed by a deconvolution or other 
method, there are two effects which can cause the estimate of the mean shape to be biased in the 
direction of the original PSF orientation. The first is a selection bias, first noted by K00. The 
second is a measurement bias which arises in the presence of anisotropic noise or anisotropic PSF, 
which we term the centroid bias. K00 discusses an effect whereby the errors in centroid appear in 
second order as biases in the ellipticity, and concludes that such errors are probably negligible. We 
show below that this bias is significant, and in fact just one of a class of noise-rectification biases 
that can occur. 



8.1. The Selection Bias 

Kaiser's selection bias operates as follows in the presence of an anisotropic PSF: if the PSF is 
elongated in the x direction (e+ > 0), then objects with intrinsic shape e+ < cover a larger area 
after PSF convolution than do objects with intrinsic e+ > 0. On the observed image, therefore, such 
objects have both lower surface brightness and lower significance v. As most detection algorithms 
involve some cut in surface brightness or u, the detected population will be biased toward e+ > 0. 
The mean e + of the population will be biased, therefore, even if all the detected objects are perfectly 
corrected for the PSF anisotropy. K00 demonstrates that the bias will scale roughly as e*o~1/o~fv 2 , 
where v is the detected significance. 

The selection bias may be defeated by careful definition of the sample of target galaxies. The 
key is to produce some significance statistic v which is independent of the shape of the object. 
For an image of flux / covering area A in an image with white noise density n, the significance is 
normally u oc f/y/nA. But the observed area A is shape-dependent in an anisotropic fashion if the 
PSF is elliptical. We may instead define v oc f/y/nA, where A is the object's area on a version of 
the image which has been deconvolved, or at least had the PSF rounded by a convolution kernel, 
as described in the previous section. We then define the target sample by some cutoff limit v m \ n to 
this "isotropic" significance. 

Detection algorithms generally have some statistic which is used as a cutoff. For KSB's imcat, 
this is v. For SExtractor and FOCAS, this is the number of pixels N above some isophotal threshold 
in a filtered version of the image. To eliminate the selection bias, we make a scatter plot of the 
detection statistic, e.g. N, vs. the "isotropic" significance v. The selection bias is eliminated if 
we choose our cutoff Dq sufficiently high that no members of the selected population approach the 
cutoff iV m i n of the detection statistic. This is illustrated in Figure 7. In this way we insure that 
our selected population is free of any anisotropic selection criterion that may have been created by 
the original detection algorithm. 
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Fig. 7. — The scatter plot of isotropic significance v vs pixels above the SExtractor thresholds, iV p i x , indicates 
the criteria necessary to avoid the selection bias. To avoid selection bias, we want the galaxy selection to be 
entirely on the basis of the shape-independent statistic v, and not influenced by the potentially shape-biased 
SExtractor threshold N m - ln = 5 [some galaxies with N < 5 are present due to object splitting]. We must 
therefore set our is m i n threshold at ss f 0, which is the lowest value for which the galaxy population does not 
extend to N < N m i n . Galaxies in the shaded region are therefore excluded. 
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8.2. Centroid Bias 

The centroid bias, discussed in K00, can be qualitatively explained as follows: an error Sx in 
centroid determination along the x axis inflates the measured second moment I xx by some amount 
oc Sx 2 , whereas a y-axis error inflates I yy . If the centroid errors are isotropic, the mean effect upon 
e + is zero. But if the x centroid errors (5x 2 ) exceed those in y, there is a net tendency to measure 
a positive e+. If the PSF has e\ > 0, then the PSF — and consequently the mean galaxy image — is 
more extended in the x direction, and hence centroids are less accurate in x than in y. The noise 
in centroid estimation therefore pushes the mean measured e in the direction of e*. The centroid 
errors scale as i/ -1 , and therefore the centroid bias in e scales roughly as e*v~ 2 . Faint galaxies at 
low v have a much larger centroid bias than the bright, high-zv stars used as PSF templates, so the 
effect is not properly removed by the KSB "smear polarizability" corrections. 

Deconvolving the image eliminates the PSF anisotropy but does not eliminate the centroid 
bias, because the deconvolved image will have anisotropic noise. There will be more noise power 
along k x than k y after deconvolution, hence there will still be anisotropic centroid errors and a bias 
toward the original e*. 

The same situation arises if we apply a convolution kernel to symmetrize the PSF. This kernel 
will smooth the image slightly in the direction perpendicular to the PSF, reducing the noise level 
in that direction somewhat. The convolved image will therefore again have a higher centroid error 
along the original PSF axis, and a consequent bias in e. 

A quantitative understanding of the centroid bias can be gained from Equation (6-53), the 
second-order expression for r\ in the simultaneous solution for shape and centroid (and possibly 
size). Let us presume that each b pq has an true value plus some measurement noise 5b pq . The 
measurement noise has mean value of zero since b pq is a linear function of the intensity. If the object 
is intrinsically round, centered, and the weight is properly sized, then we have &oo = ^io = ^20 = 0, 
and we would indeed measure rj = 0. In the presence of noise, however, the second-order terms 
in Equation (6-53) — for instance the term oc fcfo/^oo m ^ ne second line — may have non-zero mean 
values. If (frfo) = Var(&io) 7^ 0, there will be a non-zero (77). 

Note that there are a number of second-order terms in the r\ solution, arising not just from 
the solution for object center, but also from the solutions for object size and ellipticity. We apply 
the term "centroid bias" to all of these effects. 

According to Equation (6-17), Var(6io) = in the presence of a white noise spectrum P{k) = n 
(recall that 610 is a complex number). If we have a PSF of size 0+ and ellipticity e* along the x 
axis, we may attempt to round out the kernel by smoothing along the y axis a little bit. This will 
produce a noise power spectrum P'(k) ~ n(l — 2e*a 2 k 2 ). The k 2 term in the noise power produces 
a non-zero value for Var(&io), and also non-zero covariances between all the other second-order 
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elements of Equation (6-53). The leading term in the expression for centroid bias takes the form 

e = K— — « 2 ' i 8 " 1 

where K is a constant, e D is the ellipticity measured on the convolved image, a is the size observed 
on this image, and R is the resolution factor of §6. Note that the bias is in the measured shape, 
before correction for PSF dilution, which would add a factor R in the denominator. The value 
of K depends upon some higher-order moments of the typical galaxy, and upon the details of the 
PSF-correction and measurement procedure — in some cases K can be negative. Note also that the 
functional dependence of the centroid bias is essentially the same as that for the selection bias, 
hence they are difficult to distinguish. 

Figure 8 demonstrates the existence of the centroid bias in a very numerical simulation, in 
which we convolve a circular Gaussian source, with an elliptical Gaussian weight, then measure 
the ellipticity and centroid with a fixed-size circular Gaussian weight. It is clear that the mean 
measured ellipticity depends upon the significance {i.e. S/N level) as described by Equation (8- 
1) with K = —2. If the weight size is iterated to maximize significance, the bias increases to 
K ~ —6. For other measurement algorithms or galaxy shapes, K will differ, so we must determine 
K empirically. 

The bias in Equation (8-1) will be equal to the shape-noise uncertainty in the mean of N 
galaxies when 



v 1 y/N \Ke* 1 - R / 

For typical PSF ellipticities and K values, with galaxies somewhat resolved (R ~ 0.5) and detected 
at v = 10, we find the bias equal to the shape noise when N ~ 10 4 . The effect clearly cannot be 
ignored in present-day surveys. If the PSF ellipticity varies, then the bias will induce false power 
into the distortion power spectrum, with a total fluctuation power of ~ ((e*) 2 )K 2 /v 4 , for R ~ 0.5. 
If we take the RMS value of e* to be about 3%, take \K\ ~ 5, and v = 8, we see that the bias power 
is about 5% of the typical cosmic signal, for which the RMS distortion is ~ 1%. In unfavorable 
parts of the power spectrum, the ratio of PSF bias power to cosmic signal will be worse; hence the 
bias cannot be ignored for cosmic shear studies which aim to move beyond mere detection to true 
precision measurements. 

There are several possible strategies for defeating the centroid bias. The simplest is to find K 
empirically, then apply an appropriate bias correction to each object using Equation (8-1). Another 
approach, useful if one has applied a PSF-rounding kernel to the image, is to add noise back to the 
image to recreate the original isotropic noise spectrum. If both the PSF and the noise spectrum are 
isotropic, there is no way the mean shape can be biased. Unfortunately, creating the appropriate 
noise field to add to the image is computationally expensive for all but the simplest convolution 
kernels. 

Symmetrizing the noise is easier if we have used the Laguerre decomposition method to de- 
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Fig. 8. — The effect of centroid bias is demonstrated by a simple numerical test in which a circular Gaussian 
in measured in the presence of an elliptical PSF and white noise. High-significance detections are observed 
at the correct ellipticity and would yield e = when corrected for the PSF cllipticity. But the measured 
ellipticity drops as the significance v decreases. The test results are well described by Equation (8-1) with 
K = —2, which is plotted as the solid (dashed) line for the tests with e*(l — R) = 0.2 (0.1). 
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convolve the images. If we have properly propagated the original covariance matrix for the b pq , 
then we know Var(6io) — and all the other relevant covariances — in the deconvolved image. We can 
add noise to the deconvolved b pq elements in order to zero out the asymmetric elements of the 
covariance matrix. 



9. Procedures for Shape Measurement with PSF Correction 

We now have all the necessary tools for a procedure for measuring shapes and shear in the 
presence of a PSF convolution. Figure 9 is a flowchart for this procedure, and we elaborate on 
each step below. First note that there are two branch points. The first is the decision whether 
to measure the shapes on a summed image, or to measure PSF-corrected moment information 
from each exposure and combine the moments. The former is easier if all exposures cover the 
same sky area, but the latter is recommended whenever the coverage is interlaced on the sky. The 
second branch point is deciding whether to use an analytic or a kernel correction for PSF bias. 
The remainder of this section delineates the overall data reduction procedure for weak lensing 
measurements. 



9.1. Exposure Processing 

The steps listed in this section must be performed on each exposure. For mosaic detectors, 
these steps should be done separately for each channel of the detector if the PSF may change 
abruptly at the boundaries between detectors. 

1. Bias Subtraction and Flat- Fielding: These can be done in the usual fashion. Note that field 
distortion leads to significant variation in pixel area in many CCD cameras. If the flat-field 
image is, as usual, obtained by exposing a source of uniform brightness, then the flattened 
pixel data are a properly calibrated map of the surface brightness of the sky rather than a 
map of the flux from the sky. This is actually what we want, though it means that simple 
aperture photometry will lead to incorrect flux estimates. 

2. Object Detection: Available detection packages such as FOCAS (Valdes 1982), SExtractor 
(Bertin <fc Arnouts 1996), or ProFit (P. Fischer, private communication) can be used to 
identify all objects on the exposure with significance v > 3 and produce a catalog with 
preliminary position, size, and ellipticity estimation. Each of these packages also produces a 
useful estimate of total magnitude, which we will preserve. The choice of object detection 
package is not critical, because the objects that will have useful shape information have 
sufficiently high significance that any decent detection algorithm will work. 

3. Preliminary Measurement: Our shape-measurement algorithm is applied to each detected 
object in order to obtain an accurate Gaussian- weighted estimate of centroid, size (i.e. the a 
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Fig. 9. — The data flow from raw images through a final lensing distortion map or statistic is illustrated here, 
with each step described in §9. Shaded regions represent steps that may not be used in all circumstances. 
Implementation details will be presented in Jarvis et al. (2002) 
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which maximizes significance), and shape. Objects with saturated or bad pixels are flagged 
at this point. 

4. Registration: The map from pixel coordinates x p to the global celestial coordinate grid x is 
established by fitting the measured centroid of the bright objects to a collective catalog (e.g. 
the average of all individual exposures' catalogs) and/or an astrometric catalog. This step 
is critical as even slight misregistrations will produce a systematic coherent ellipticity to the 
summed images. The photometric offsets between images are determined at this point as 
well. 

5. Image Summation (optional): The exposures can now be remapped to a common grid and 
photometric scale and summed to give a deeper image. If the flattened images are surface- 
brightness maps as described in Step 1, then a simple interpolation can be done when remap- 
ping the images — no Jacobian flux corrections are needed. Recall the caveats about image 
combination in §4.1. If shapes are to be measured from the summed image, then we must 
at this point reconstruct the preliminary catalog by executing Steps 2 and 3 on the summed 
image. 

6. Identification of Stars: As many unsaturated stellar images as possible should be identified on 
the image for use in PSF determination. Stars are typically identified on the size-magnitude 
plane. This is a difficult task to execute without excessive human intervention; we want our 
star-finding algorithm to be flexible enough to identify stars successfully in the presence of 
a spatially varying PSF, but strict enough to insure that galaxies are not counted as stars. 
Failure on either count will lead to an erroneous estimate of the PSF size and shape on some 
part of the exposure, which will lead to a false feature in the shear field. Our algorithms for 
star identification are described in Jarvis (2002). The algorithms are particularly suited to 
identifying stars in the presence of a spatially varying PSF size. 

7. PSF Measurement: The shape-measurement algorithms are now applied to the identified 
stars. The integrals in Equation (6-16) give b* at the location of each star. Note that the 
integral can readily be executed in the global coordinate system x because we know the map 
from pixel coordinates, and because pixel data are already in terms of surface brightness /. 
Thus the effects of image distortion are removed at this step. 

8. PSF Interpolation: With the PSF b* measured at the locations of the stars, we need to fit a 
function b*(x) which describes the PSF at any location on the image. We use a polynomial 
to describe the variation of each 6* across the image. Note that for complex or undersampled 
PSFs (such as those on HST), interpolation of the b* components is much easier and more 
accurate than interpolating a pixel-wise representation. In any case this step is again critical: 
the interpolation scheme must be flexible enough to follow PSF variation, but must remain 
well-behaved at the image edges and in regions where PSF stars are sparse. There must 
also be some form of outlier rejection for PSF stars that have cosmic rays or near neighbors 
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contaminating the measurement, but a true excursion of the PSF in some part of the image 
must not be rejected. 

9. PSF-Rounding Kernel (optional): At this point we derive and apply the kernel to remove the 
anisotropy from the PSF, if desired. Our method for doing so is described in §7. After the 
convolution it is necessary to repeat Steps 7 and 8 so as to have an updated map of the PSF. 



9.2. Object Measurement 

We next outline the procedures which yield a shape estimate for every object in the field. 

10. Creation of Target Catalog: A list of all targets for shape measurement must be compiled. 
If there are < 5 exposures of each object, then a master target list can be produced by 
taking the union of the individual exposure catalogs produced in Step 2. One can use a 
very low detection threshold on the individual exposures, then guard against noise detections 
by demanding coincidence on 2 or more exposures. If there are > 5 exposures per object, 
then it will be necessary to run a detection algorithm on a summed image in order to find 
all the potentially useful target galaxies. The target catalog should include an estimate of 
the centroid (in global coordinates) and some observed size s of each object. To avoid 
Kaiser's selection bias (§8), the criterion for acceptance into the catalog must be some shape- 
independent statistic. 

The following steps are performed for each object in the target catalog: 

11. Measurement of Observed Moments: Given the global coordinates of the object, one can de- 
termine all of the exposures on which it should appear. Note that for low-significance objects, 
this may include exposures for which the object was missed in the preliminary detection of 
Step 2. We then use Equation (6-16) to measure the Laguerre expansion b° of the image as 
observed on each exposure. Once again the integrals are performed in the global coordinate 
system to remove the effects of optical distortions. The integrals require choice of a centroid 
and a size parameter a for the basis functions ip". We may use the approximate centroid 
from the target catalog, and set the weight size a equal to the typical observed object size s 
from the target catalog. The covariance matrix for the b° vector has the diagonal form given 
in Equation (6-17) for sky-dominated galaxies. Measurements contaminated by saturated 
pixels or other defects may be rejected at this step. 

12. Correction for PSF Effects: There are two methods available here. The first sequence may 
be used if the rounding kernel has been applied to remove the PSF bias. In this case, if the 
significance on the individual exposure is > 3, we may proceed as follows: 
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(a) Shape Measurement: The image can be shifted, dilated, and sheared until it passes the 
centroid, significance, and roundness criteria embodied by Equations (6-19)-(6-21). This 
yields an optimal measure of the observed ellipticity e°. If the object is of significance 
u o It 3 on the individual exposure, then we can solve for centroid and size on each 
exposure independently and still ignore the higher-order terms in the uncertainty of 
Equation (4-1). 

(b) Centroid Bias Correction: Each measured ellipticity must be corrected for the centroid 
bias (§8) using an empirical value of K in Equation (8-1). 

(c) Dilution Correction: The PSF is interpolated to the position of this object to determine 
the resolution parameter R (Equation [6-1]) or the higher-order version derived in Ap- 
pendix C. The correction for dilution to give the image-plane (pre-seeing) ellipticity is 
then simply e* = e°/R. The uncertainty ellipse for e (Equation [3-30]) is therefore also 
scaled by 1/R. Note that it is possible to obtain e l > 1 if the noise makes the target 
appear smaller than the PSF. 

(d) Averaging of Exposures: The e* from the collection of exposures are averaged using 
the weighting procedures of §4.2. Some form of outlier rejection is necessary to remove 
objects contaminated by cosmic rays or other defects. 

If we are to perform the PSF bias correction analytically, or if the significance per exposure 
is low, then it is best to average the deconvolved moments rather than deriving e* for each 
exposure: 

(a) Deconvolution of Moments: The PSF decomposition b* is interpolated to the position 
of the object, and the form of the convolution matrix C(b*) is calculated using the 
coefficients in Appendix B. This matrix is truncated at some order p + q < N and 
inverted; the deconvolved (pre-seeing) Laguerre decomposition is b* = C _1 • b°. Since 
this is a linear operation on b, the covariance matrix for b* can be propagated from the 
simple diagonal covariance matrix for b° (Equation [6-17]). There is a subtlety involved 
in the choice of weight scale er : in §6.2 we determined that the ideal weight scale for 
the deconvolved image is erf = sf + o\ = s^, where st and s are the sizes of the pre- and 
post-seeing objects, and cr* is the size of the PSF. The typical value of s Q was placed in 
the target catalog in Step 10. According to Equation (6-63), our deconvolution formulae 
are simplest if we choose the weight scale by 

2 2,22,2 / n 1 \ 

a =o- i +0-+ = s + cr*. (9-1) 

Thus in Step 11 we in fact want to use a weight scale somewhat larger than the s Q which 
maximizes the significance of detection. 

(b) Combination of Moments: From each exposure we have estimated a b* deconvolved 
moment vector, with known covariance matrix. We average these vectors to obtain a 
single best estimate of b\ We have a choice of weights to apply in producing the average; 
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an obvious choice is to weight each exposure inversely with Cov(&2o&2o)> shrce &20 carries 
most of the shape information. Again some form of outlier rejection is necessary at this 
step. 

(c) Symmetrization of Noise: Certain elements of the propagated covariance matrix for the 
b vector must be zero in order to avoid noise anisotropies that produce centroid bias 
(§8). This can be achieved by selectively adding Gaussian noise to various elements of 
the combined moment vector. 

(d) Determination of Shape: Given the average of b* from all exposures, we find the trans- 
lation and shear that must be applied to satisfy the centroid and roundness conditions 
bio = ^20 = 0. The formulae of Appendix A are used for this. The covariance matrix for 
b* can be propagated through the transformations as well. The uncertainty in the shape 
rj is then the square root of Var(620^2o)/(^oo)- Note that we do not want to maximize 
the significance by dilating to set bu = 0. Our optimization criterion is that the shape 
have minimal error after our transformations. One could dilate the image to satisfy this 
desire, but our choice of cij should already have us close to the optimum, according to 
the arguments of §6.2. 

After either of these procedures, we have a measure of the deconvolved shape along with its 
uncertainty. 

13. Combination of Different Wavelengths: If we have imaged the field in a variety of filters, 
then we will have obtained a shape measurement in each filter. The galaxy's shape and 
moments may depend upon wavelength so we don't want to average together images or 
moments measured in different filters. We can, however, use the methods of §4.2 to produce 
a wavelength-averaged e l which is maximally sensitive to shear. Weighting each filter by 
the error in its measured shape insures that we obtain the most sensitivity from each galaxy 
regardless of its color. 

After completion of all these steps, we have a catalog of all objects in the field, specifying their 
location, magnitude, optimal shape measurement, and shape uncertainty. 



9.3. Determination of Shear 

With the shape catalog in hand we are close to the scientific goals. The remaining step is 

14. Generation of Shear Data: The target galaxies are binned by position, etc., into subsets for 
which we wish to determine a shear. The shapes (and their uncertainties) may have to be 
rotated, e.g. into tangential coordinates about some mass center, depending upon the shear 
statistic under study. The formulae of §5 take the collection of shapes and uncertainties and 
allow us to create an optimal shear estimate, as well as to propagate the uncertainties in 
shape to the shear measurement. 
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10. Conclusions 

We have attempted to produce, as rigorously as possible, an end-to-end methodology for 
measuring gravitational distortion that has optimally low noise, with calibration and noise levels 
derivable entirely from the observations themselves. We have succeeded in many, but not all, 
aspects of the problem: 

• The measurement of individual galaxy shapes appears to be optimal and has traceable noise 
characteristics; indeed there even appears to be a straightforward way to handle undersamplcd 
data and retain proper covariance information for the intrinsic Laguerre coefficients b l pq (§6.4). 
The Gaussian weights underlying the Laguerre decomposition are nearly optimal for sky- 
dominated exponential-profile galaxies. We may wish to reexamine this scheme for the case 
of when the galaxy is brighter than the sky background. 

• The correction of measured moments for the distorting and diluting effects of the PSF can 
be effected to arbitrary accuracy using the Laguerre-decomposition methods. There will 
be a tradeoff between elimination of systematic errors — which pushes toward inclusion of 
higher-order terms in the deconvolution — and the minimization of measurement noise from 
high-frequency terms. It is not clear if the Laguerre method is optimal with, for example, 
Airy PSFs with sharp cutoffs in k space; but the method should be better than those yet 
applied. 

• We have identified methods to work around two sources of bias that arise from PSF ellipticities 
even in the presence of perfect deconvolution. 

• Measurements of galaxy shapes from different exposures or filter bands can be optimally 
combined with standard least-squares techniques, since we know the e uncertainty from each 
individual exposure. 

• The determination of lensing distortion from the ensemble of galaxy shapes has a straight- 
forward, optimal, calibratable solution in the case of no measurement noise. In the presence 
of measurement noise, however, it is necessary to have some knowledge of the noiseless shape 
distribution to get the calibration factor exactly correct. Our approximations, however, seem 
to suffice to get accuracies of 5-10%. 

A detailed performance comparison of our methods with other authors' is beyond the scope of 
the paper. We can guess, however, that the reduction of measurement noise relative to a carefully 
weighted implementation of KSB will be slight, perhaps a factor of 1.5. But our methods, like 
those of K00 and Refregier (2001), are formally valid for any reasonable PSF, and hence we expect 
to have much-reduced systematic errors. Indeed it is likely that the biases of §8 have not been 
extensively tested because they were lost under the larger errors in PSF correction (as also noted 
by K00). 
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Our methods share elements with many previous proposals. An aesthetic difference is that we 
retain the geometric meaning of ellipticity by, in effect, using adaptively shaped weights instead 
of fixed circular weights. This makes the "shear polarizability" a purely geometric effect. As a 
consequence we can examine the P{e) distributions and find ways to exploit surface brightness, 622, 
or color information to separate the spheroid galaxies and weight them more heavily to reduce the 
shape noise. 

It will be of interest to see how our method compares to the commutator method of K00. One 
would hope that the two independent methods could be applied to the same dataset and yield the 
same results, bolstering our confidence in these very difficult measurements. 

In a succeeding paper (Jarvis et al. 2002) we will present some of the implementation details 
for the analytical methods here, and test the methods on real and simulated data. Fischer et 
al. (2000); Smith et al. (2000); Wittman et al. (2000) make use of portions of this methodology, 
so the systematic-error tests and calibration tests in those papers already serve as demonstra- 
tions. Upcoming precision measurements of cosmic shear will make even greater demands upon the 
systematic-error reduction and accurate calibration that our methods offer. 

This work was supported by grant AST-9624592 from the National Science Foundation. We 
would like to thank Phil Fischer, Deano Smith, Tony Tyson, Jordi Miralda-Escude, and Dave 
Wittman for many discussions on these methodologies, and help in implemented various image- 
processing algorithms. Thanks also to all our collaborators who have waited several years for a 
coherent explanation of the methods we are using. 
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A. Point Transformation for Laguerre Expansions 

In §6.3.2 we derived the mapping matrices for the vector b of Laguerre-expansion coefficients 
when the underlying image / is transformed by an infinitesimal translation, dilation, or shear. In 
this Appendix we derive the coefficients of the mapping matrices for finite transformations. We 
have implemented these transformations as methods for C++ classes that represent the Laguerre 
expansions. 



A.l. Translation 

We defined the transformation T 2 on the image by 

T z I(x,y) = I(x -x ,y -y ) (Al) 

z ee f°±M. (A2) 

a 

We first define a translated raising operator ap by 

ap(T z I) = T z {apl) (A3) 

=> ap = <-|. (A4) 

The second line is apparent from examining the form of the raising operators in Equation (6-22). 
We decompose the images into our eigenfunctions as I = J2 b pq ipp q and T Z I = Y^ b' pq ipp q ; then we 
can express T z as a matrix operation on b: 

(A5) 
(A6) 

(A7) 
(A8) 



(A9) 

Applying the translated raising operator (A4) to the definition (A7) of T v ? i yields the recursion 
relation 



b' = 


T 2 b 


b p , q , = 


Z^ p'q'"P1 


TW& = 


}2 T p^p'q'- 


1 p'q' - 


a 2 fd 2 x{T z r pq Wp lql . 


i the first coefficient 


T-.00 


' = e"l z l 2 / 4 . 



P' %U )q > - i 2JJ = y/p + lTfi 1 ". (MO) 



2 
The same procedure using the q raising operator and the lowering operators gives the recursions 



^t- i)-^ P v = v^ttt;/; +1 > (ah) 



p'(q'-l) 2 P'l' V P'l' 

V¥TlTff +1)ql - 1^, = ^T^ 1)q (A12) 

Z rppq _ f^rrPiq- 1 ) 



V¥+ 1T ?(q' + l)-0 T p'q' = ^ T Tq> ^ 
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These two relations allow us to generate any T^?, recursively from Tqq . In particular, using the 
first 2 we have 

T m = (-i/y^ p-g^iv^ (A14) 



00 VpW- 

One can derive a closed-form expression for the general T^?, using further recursion, but the ex- 
pression involves a double sum and is not particularly illuminating. The most efficient algorithm 
for computing all the coefficients is to note that the generator for translation (Equation [6-38]) sep- 
arates into p- and ^-dependent components. The general matrix element can therefore be expressed 
as 

T™, = f(p,p')f(q,q'), (A15) 

with the above recursion relations for TZ q , leading to 

/(P,0) = ^! e -MV8 ; (A16) 



f(p,P+l) = [Vpf(p-hp') + ^f(p,P > )\/VP T +^ (A17) 



A.2. Dilation 

The dilation operation is defined by 

D M /(x,y)=/(e-^,e-' i y) (A18) 

As for the translation, we can define a dilated raising operator and use it to derive a recursion 
relation for the coefficients of D: 

o*t(D A1 /) = B^apl) (A19) 

=4> aZ * = cosh /j, ap — sinh /j, a q (A20) 

,(p+i)g 



=> V / P+l^ ; * = cosh M ^£»5_ 1)9 , - S inh/iV7 + Ti^ (g , +1) (A21) 

Using the recursion operator and its q equivalent, we can generate any desired coefficient from D®®. 
With the direct integration analogous to Equation (A8), we derive 

Df q = e^sech^tanh nf S pq . (A22) 

In deriving this we make use of the identity (Abramowitz & Stegun 1965) 

LfO(as) = £ ( 9 ^V U " ar fc 4 m) (x). (A23) 

fe=0 ^^ ' 

In fact with this identity one can derive any D v 3 , by direct integration, but the closed form is a 
sum over k that is again not particularly useful, as the recursion relation (A21) is a faster way to 
generate the coefficients. 
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A.3. Shear 

A shear r\ oriented on the x axis gives the transformation 

S v I(x,y)=I(e- r '/ 2 x,e^ 2 y). (A24) 

When we define the transformed raising operator dp as for the translation and dilation, we find 
the operator and consequent recursion relation to be 

=► ap = cosh(r//2) a^ - sinh(r ? /2) a a p (A25) 

=* v^TT^ 9 = cosh(7//2)^SW_ 1)9/ -siiili(7//2)v^+lS^ +w . (A26) 

The recursion again simplifies by noting that the shear generator in Equation (6-33) separates into p 
operands and q operands. Hence the matrix elements must be expressible as S^f , = f(p,p')f(q, </)• 
We can determine the function f(p,0) by direct integration of the S p matrix element, which yields 
(for p even) 

QP o VP 1 - w /0 ^ -tanh(r//2) y /2 
°° = fo/2)T ^ ^ 1 2 J ( ^ 



U 



/p!g! 



tanh(^/2)V +,)/2 



s « = w» sech(,,/2 H^^ ' (A28> 



The coefficients vanish if p or q are odd. The recursion relation then generates any desired coeffi- 
cient. For a shear oriented at position angle f3, the coefficients acquire an additional phase factor 
exp[i(p' - q' -p + q)0\. 



B. Convolution of Laguerre Expansions 

We wish to derive the coefficients which express the convolution of two eigenfunctions as a 
new sum over eigenfunctions, as defined in Equation (6-57). This will be easier if we work in k 
space, with the convolution turned into a multiplication as in Equation (6-59), with the fc-space 
eigenfunctions given in Equation (6-62). 

A rapidly computable recursive formulation of the coefficients is again derivable from the 
raising operator. From the form of the fc-space raising operator (6-61) we see that if o 2 = o 2 + a 2 , 
then 

(Jocf/lo = {aiafli^ + Iiia^Q (Bl) 

=* ^v^+lC^*^ 1 )** = 'oyfcC^fc-aiy/^Cto* 1 )****. (B2) 

An equivalent manipulation with the lowering operator yields the recursion relation 

n frT n PiqiP * q * — rr /7T r ( ~P^ l ^ qiP * q * -I- a /7T r Piq ^ p *^ q * ('R'W 
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These recursion relations, and their q equivalent, will allow us to derive any coefficient if we know 
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Volo 1 



i.e. if we can calculate the effect of multiplication by the Gaussian ip™. This is straight 



forward if we recall that m = rrii for m± = 0, and use Equation (A23): 
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erf, 



(B4) 

(B5) 
(B6) 



The parameter D is the "deconvolution ratio," with D = 1 being the limit of perfect resolution and 
D = meaning a PSF much broader than the image. This expression and the recursion relations 
are both separable into p- and (/-dependent terms, so we can simplify the computation of the matrix 
elements by using the expression 



nPiqiP*q* 
p q 



A 



2^ 



PolAl 



G(p ,pi,p* 



qolAl 



G(q ,qi,q*) 



Pi + p* - Po = qi + q* - q > 0. 



(B7) 
(B8) 



Terms for which the conditions in the second line are not met are zero. The recursion relations and 
specific values given above for C can be recast for the function G as follows: 



G(p + l,Pi,p± 



—G(p ,Pi - l,p*) + —G(p ,p i ,p i , - 1) 



.UP* 



Pi+P-> 
Pi 



p* 



(B9) 
(BIO) 



The symmetry between initial and PSF images is clear in these equations. There is a consequent 
closed form for G(p ,pi,p+), but it is again not particularly illuminating, and the recursive form is 
stable and faster for computation. 

In the case where the PSF is a unit-flux Gaussian, Equation (B5) can be used to give the 
observed b° in terms of the intrinsic decomposition: 



b°p q =E D{P0+q ° )/2 < 
3=0 



Po + j\ (q + j 

Po 



( l - D )% +i)iq +i) (Po>Vo). 



(Bll) 



Note that the convolution matrix C is in this case block-diagonal, as the m states do not mix, and 
also upper triangular, as p ,q < Pi,qi for non-zero elements. The inverse (deconvolution) matrix 
can in this case be expressed in closed form: 



b 



p,q, 



3=0 



Pi + j\ (Qi+j 

Pi J\ qi 



1- D 
D 



(-l)'fr' 



(p»+j')(?»+j)" 



(B12) 
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Second-order Formulae for PSF Dilution 



Here we describe a refinement to the resolution parameter R defined in Equation (6-1). The 
observed ellipticity is R times the true ellipticity in the special cases of unweighted moments or 
when the PSF is a circular Gaussian and the galaxy an elliptical Gaussian. In these cases the factor 
R can be expressed as 



s 2 



R = 1-H*, (CI) 

s 2 = {(x 2 + y 2 )/2) = a 2 cosh??. (C2) 



In the second line we have assumed that the methods of §3 have been used to shear the object by 
r] to produce something that is "round" under a Gaussian weight of optimal size a. 

Here we derive a form for R which is applicable to the case in which the galaxy has homologous 
elliptical isophotes with i|Cl, but with an arbitrary radial profile. In the language of our Laguerre 
coefficients, we have b\\ = by proper choice of a, and b pp are arbitrary for p > 2. By Equation (6- 
34) this slightly out-of-round object has 620 = v2?7(&oo — &22)/4. All odd-indexed coefficients arc 
zero. 

We also take the PSF size <r* to be small compared to the intrinsic object size a*. The action 
of this small isotropic convolution C CTvt on the moments is the same as a transformation in which 
the image I is replaced by the average of two versions displaced by ±cr* in the x direction, followed 
by a similar infinitesimal spread in the y direction. Defining z = (T±/a and using a second-order 
version of the generator Equation (6-38), we can show that 

C CTi = i (T z + T_ 2 ) (T iz + T_ i2 ) (C3) 



z 2 

1 + T 



a P a\ + a p a q - p - q - 1 



(C4) 



When this convolution acts upon our original slightly-elliptical object, the resulting object has 
Laguerre coefficients b' with 

b' 00 = b 00 (l-z 2 /2) (C5) 

b'20 = b 20 (l-3z 2 /2) (C6) 

b' 22 = b 22 (l-5z 2 /2) (C7) 

We then need the value rf of the shear that will make this new object appear round. According to 
Equation (6-47), this will give 

, r?(6oo-6 22 )(l-3z 2 /2) 

V 1 - z 2 /2 - 622(1 - 5^2/2) [ ' 



R^^ = l-z 2 (l + - u \ +0(z 4 ) (C9) 

ooo — b 22 _ 
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oi 



(CIO) 



g2 s ^60^622 (cn) 

»00 + »22 

Note that the kurtosis measure 04 defined in Equation (3-26) is the same as 622 /&00 which appears 
here. We make the ansatz that the correct form for R in the case of finite dilution or finite e is the 
Gaussian form Equation (CI) with the kurtosis term added: 

ei = e /R, (C12) 

R = 1-%, (C13) 

s 2 = i^VcoshTj. (C14) 

1 + a 4 

Note that we apply the kurtosis correction to the PSF size measure s 2 in the same way as for the 
object, to give a well-behaved correction for poorly-resolved objects. 
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